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Abstract 

Small  unmanned  aerial  systems  are  being  designed  to  emulate  the  flapping  kinematics 
of  insects  and  birds  which  show  superior  control  in  slow  speed  regimes  compared  to 
fixed  wing  or  rotorcraft  aircraft.  The  flight  of  flapping  wing  vehicles  is  characterized  by 
aeroelastic  effects.  Wing  membranes  are  highly  flexible  with  large  deformations  during 
the  stroke.  Predicting  the  aeroelastic  response  of  a  flapping  wing  is  an  essential  step 
towards  predicting  performance  characteristics  of  different  designs.  Most  research  has 
been  dedicated  towards  understanding  the  aerodynamic  side  of  the  aeroelastic  response 
with  minimal  effort  spent  towards  validating  the  structural  response.  A  finite  element 
model  of  a  wing  from  a  commercial  flapping  wing  vehicle  was  created  to  validate  the 
structural  response.  Vacuum  testing  allowed  the  isolation  of  the  inertial  response  for  a 
direct  comparison  to  the  finite  element  model.  Model  tuning  was  accomplished  through 
modal  testing.  The  finite  element  model  was  run  through  a  nonlinear  dynamic  simulation 
to  match  the  test  article  flapping  motion.  Wing  tip  displacement  amplitude  was  matched 
to  within  8%.  The  membrane  kinematics  of  the  finite  element  model  were  similar  to  the 
vacuum  test  article  but  overall  membrane  deflections  predicted  by  the  finite  element 
solver  were  less  than  observed  deflections  seen  in  the  vacuum.  Modeling  the  wing 
membrane  dynamics  proved  to  be  difficult  because  of  the  nonlinearities  associated  with 
the  flexible  membrane.  This  research  shows  that  significant  focus  must  be  placed  on 
validating  the  structural  side  of  a  flexible  structure  in  order  to  correctly  model  the 
complete  aeroelastic  response. 
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FINITE  ELEMENT  ANALYSIS  OF  A  HIGHLY  FLEXIBLE  FLAPPING  WING 


I.  Introduction 


1.1  Why  Small  Unmanned  Aerial  Systems 

When  manned  flight  first  became  possible,  researchers  focused  on  ways  to  travel 
higher,  faster,  and  with  larger  payloads.  Now  that  technology  has  allowed  the 
miniaturization  of  payloads,  research  agencies  are  shifting  their  focus  to  develop  smaller, 
slower  aircraft  to  accomplish  complicated  missions.  In  1996  the  Defense  Advance 
Research  Projects  Agency  (DARPA)  established  goals  of  a  Small  Unmanned  Aerial 
System  (SUAS)  program  listed  in  Table  1  [1].  The  Air  Force  Research  Laboratory 
(AFRL)  set  out  to  build  an  operational  bird  size  SUAS  by  2015  and  an  insect  size  SUAS 
by  2030.  The  problem  is  this  reduction  in  size  creates  new  challenges  for  aircraft 
designers. 


Table  1:  SUAS  Goals  set  by  DARPA  [1] 


Parameter 

SUAS  Value 

Size 

<15  cm 

Weight 

10-100  grams 

Useful  Payload 

1-18  grams 

Endurance 

20-60  minutes 

Air  Speed 

30-65  km/hr  cruise,  hover  is  tradeoff  with  endurance/range 

Range 

1-10  km 

Since  1996,  the  Department  of  Defense  (DoD)  goals  for  SUASs  have  evolved  with 
the  progress  of  technology  and  feasibility  but  the  same  design  challenges  remain.  Most 
current  operational  SUASs  are  based  on  conventional  aircraft  design.  Conventional 
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aircraft  generate  lift  by  speeding  air  flow  over  rigid  airfoil  shaped  surfaces.  To  achieve 
the  required  air  speed  over  the  lifting  surfaces,  conventional  aircraft  must  physically 
move  forward  through  the  air  with  the  aid  of  a  propulsion  system.  As  such,  conventional 
aircraft  operate  in  the  high  Reynolds  number  region  and  are  generally  unable  to  perform 
well  in  slow  operating  regions  within  a  confined  space.  The  difference  in  Reynolds 
number  regimes  is  graphically  displayed  in  Figure  1  with  the  weight  of  the  vehicle. 
Miniaturization  of  conventional  aircraft  designs  to  meet  the  future  requirements  of 
agencies  like  the  DoD  and  AFRL  is  not  feasible  because  of  the  low  Reynolds  numbers 
and  unsteady  flow  involved  in  the  envisioned  operating  region  [1,2].  Figure  2  shows 
select  current  operational  Unmanned  Aerial  Vehicles  (UAVs)  sizes  compared  to  their 
payload;  the  projected  SUAS  region  has  yet  to  be  broached  by  operation  aircraft  [3],  To 
solve  these  design  challenges,  research  has  focused  on  the  flapping  motion  of  natural 
flyers  (birds  and  insects)  to  discover  features  that  allow  them  control  in  low  Reynolds 
numbers  and  unsteady  flow  regimes  with  the  ultimate  goal  of  incorporating  those  features 
into  future  systems. 
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Figure  1.  Reynolds  number  range  for  flight  vehicles  [4]. 
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Figure  2:  Payload  versus  wingspan  for  existing  UAVs  and  proposed  SUASs.  Predicted 
capabilities  fall  within  the  trends  extrapolated  from  larger  vehicles  [3]. 
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Flapping  wing  flight  offers  many  advantages  over  conventional  aircraft  for  small, 
slow  traveling  vehicles.  As  stated  earlier,  conventional  aircraft  must  physically  travel 
forward  through  the  air  to  generate  lift  and  as  such  conventional  aircraft  can  never  hover. 
For  flapping  wing  vehicles,  the  ability  to  hover  offers  a  huge  operational  advantage  for  a 
myriad  of  different  tasks.  Flapping  wings  also  have  an  aerodynamic  advantage  over 
conventional  design.  Birds  and  insects  demonstrate  an  amazing  ability  to  maneuver 
through  the  air,  and  although  this  feat  is  not  obtainable  by  current  flapping  wing  test 
vehicles,  the  potential  maneuverability  is  far  beyond  the  best  a  conventional  aircraft 
could  ever  achieve  at  slow  speeds. 

Rotorcraft  are  another  category  of  aircraft  that  could  potentially  fill  this  small,  slow 
vehicle  requirement.  Rotorcraft  offer  the  ability  to  hover  and  can  achieve  decent 
aerodynamic  perfonnance  characteristics.  However,  rotorcraft  are  loud  and  would  not  be 
an  ideal  choice  where  stealth  is  a  must  like  surveillance  and  reconnaissance.  Even  when 
not  perfonning  stealth  missions,  the  noise  pollution  generated  from  a  rotorcraft  can  be 
damaging  to  a  local  community.  Flapping  wing  offers  two  advantages  here  in  the  stealth 
category:  they  can  be  quiet  and  inconspicuous.  Flapping  wings  vehicles  can  be  very 
quiet  especially  in  nature.  Some  insects,  like  butterflies,  cannot  be  heard  unless  they  are 
quite  literally  on  top  of  a  person.  Because  flapping  wing  vehicles  mimic  natural  flyers, 
they  have  the  potential  to  blend  in  with  the  environment.  No  astute  person  will  mistake  a 
rotorcraft  for  a  bird  or  insect.  All  the  advantages  flapping  wing  vehicles  exhibit  over 
conventional  aircraft  and  rotorcraft  makes  them  an  ideal  platform  for  future  SUAS. 
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The  flapping  wing  motion  of  natural  flyers  is  a  recent  area  of  extensive  study  because 
of  its  potential  application  to  SUASs  and  has  given  rise  to  the  production  of  flapping 
wing  small  unmanned  aerial  systems.  Before  the  mid  1990s,  virtually  no  serious  research 
was  conducted  on  the  flapping  wing  aspect  of  flight.  Therefore,  replicating  a  natural 
flyer’s  wing  motion  and  dynamic  response  is  an  area  of  on-going  practice  [5,  1].  Since 
flapping  wing  effects  are  not  fully  understood,  researchers,  including  biologists  and 
engineers  alike,  turn  to  nature  for  the  answers.  While  biologist  have  extensively  studied 
the  flying  mechanics  behind  birds  and  insects  for  quite  some  time,  the  reasons  why  these 
mechanics  work  are  not  yet  fully  known  [4]. 

Natural  flyers  such  as  birds  and  insects  operate  in  the  unsteady  low  Reynolds  number 
region  and  do  so  with  precise  control  and  agility.  By  studying  these  natural  flyers, 
researchers  aim  to  understand  what  allows  birds  and  insects  to  operate  in  low  Reynolds 
numbers  regimes  successfully  and  then  apply  those  characteristics  to  SUASs.  While 
birds  and  insects  both  use  flapping  to  maneuver  through  the  air,  birds  have  complex 
wings  with  muscles  and  joints  which  allow  them  to  actively  change  the  shape  of  their 
wing  throughout  different  flight  regimes.  Insect  wings  do  not  have  muscles  and  joints  so 
the  shape  of  their  wing  in  flight  is  caused  by  a  combination  of  inertial  and  aerodynamic 
loads  acting  on  the  surface  [5].  This  study  will  also  focus  on  passively  controlled  wings 
(like  insects)  that  contain  no  mechanism  to  actively  control  the  wing  shape  in  flight. 
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1.2  Research  Focus 


Why  is  modeling  the  motion  of  a  flapping  wing  important  to  aircraft  designers  and 
future  flapping  wing  SUASs?  To  answer  this  question,  it  is  important  to  know  an 
overview  of  how  operational  aircraft  are  designed.  Somewhere  in  the  field,  someone  or 
some  agency  has  a  mission  that  needs  to  be  accomplished.  Mission  requirements  crucial 
to  mission  success  are  formulated  based  upon  these  needs.  Mission  requirements  can  be 
very  specific  such  as  detailing  a  range  and  payload  capacity  or  they  can  be  broad 
description  of  a  task  to  be  accomplished.  System  engineers  take  these  mission 
requirements  and  translate  them  into  perfonnance  characteristics.  Performance 
characteristics  are  usually  more  quantitative  measures  like  aircraft  weight,  max  speed, 
etc.  that  have  meaning  to  an  aircraft  designer.  Conventional  aircraft  designers  know  how 
to  incorporate  design  features  into  an  aircraft  to  meet  the  performance  characteristics. 
Here,  design  features  refer  to  specific  details  of  the  aircraft  like  wingspan,  aspect  ratio, 
camber,  and  a  myriad  of  other  quantities.  The  problem  with  flapping  wing  SUASs  is  that 
it  is  not  understood  how  design  features  affect  the  performance  of  the  vehicle.  Without 
this  information  designers  cannot  build  optimal  designs  to  meet  specific  performance 
characteristics.  SUASs  can  be  successfully  built  and  flown  through  trial  and  error  even 
though  why  they  fly  is  not  understood.  Current  research  is  aimed  at  figuring  out  how 
design  features  translate  into  performance  characteristics.  By  modeling  a  flapping  wing, 
researchers  will  be  able  to  predict  performance  characteristics  for  any  wing  design,  and 
aircraft  designers  will  be  able  to  truly  optimize  designs  for  operational  use. 

Simply  put,  the  focus  of  this  research  was  to  successfully  model  the  deformations  of  a 
flapping  small  unmanned  aerial  system  wing  with  a  finite  element  (FE)  analysis.  Solely 
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using  a  FE  solver  would  greatly  simplify  the  modeling  process;  however,  FE  only 
predicts  the  structural  deformations.  Computational  Fluid  Dynamics  (CFD)  is  a  very 
powerful  but  computationally  intensive  way  to  model  the  aerodynamic  effects  of  a 
flapping  wing.  While  other  research  has  successfully  used  CFD  and  FE  analysis  to 
model  a  flapping  wing,  the  application  is  still  limited  due  to  the  complexity  of  the 
simulation  and  computational  effort  [6].  Most  SUAS  researched  is  focused  on  studying 
the  flapping  wings  of  insects  such  as  the  dragonfly,  hawkmoth,  butterfly,  etc.  but  this 
research  will  focus  on  a  commercially  built  SUAS.  Most  research  shows  that  insect 
wings  behave  aeroelastically,  but  their  wings  are  still  fairly  stiff.  This  research  will  focus 
on  a  wing  with  a  stiff  leading  edge  beam  attached  to  a  highly  flexible  membrane.  As 
such,  the  relative  deformations  observed  will  be  greater  than  previous  studies  of 
biological  inspired  SUASs. 

Modeling  the  aeroelastic  response  of  a  flexible  flapping  wing  is  very  complicated  so 
the  FE  model  was  first  verified  by  testing  in  a  vacuum  chamber.  In  the  absence  of  air,  the 
deformations  observed  result  from  the  inertial  forces  of  the  wing  and  are  directly  related 
to  the  structure  and  material  properties.  The  deformations  recorded  in  the  vacuum  were 
compared  to  the  deformations  recorded  in  the  air.  The  only  difference  between  the  two 
setups  was  the  presence  or  lack  of  aerodynamic  forces.  This  difference  allowed  for  an 
assessment  of  the  dominance  of  inertial  versus  aeroelastic  forces. 

1.3  Modeling  Efforts 

Most  modeling  efforts  have  attempted  to  validate  biological  flapping  wings  of  various 
insects.  Many  researchers  focus  solely  on  the  aerodynamic  side  of  the  aeroelastic 
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response  and  little  effort  has  gone  into  validating  the  structural  models.  The  research 
efforts  that  have  focused  on  structural  modeling  have  either  used  static  analysis  efforts  or 
analyzed  relatively  stiff  wings.  This  research  will  focus  on  the  wings  of  a  test  article 
characterized  by  billow  during  the  flap  cycle.  All  the  methods  proposed  for  this  research 
have  been  used  in  previous  studies  but  none  have  been  put  together  to  model  the  dynamic 
deformations  of  a  flexible  flapping  wing  in  this  manner. 

1.4  Chapter  Overview 

Chapter  II,  the  literary  review  section,  will  cover  the  background  behind  flapping 
wing  motion  and  review  many  research  efforts  to  model  the  structural  and  aerodynamic 
response.  Chapter  III  details  the  test  article  used  for  testing,  laboratory  setup,  the 
methodology  used  in  the  research,  and  the  structural  model.  Chapter  IV  provides  the 
results  and  analysis  of  all  tests.  Chapter  V  draws  the  final  conclusions  on  the  research 
and  provides  recommendations  for  future  work. 

1.5  Preview  of  the  Results 

Correctly  modeling  all  the  intricacies  of  flapping  wing  motion  proved  to  be  a  difficult 
task.  The  dynamic  nature  of  the  flapping  motion  coupled  with  the  nonlinearity  of  the 
wing  membrane  meant  precise  modeling  of  the  problem  was  required.  The  FE  wing  tip 
displacement  amplitude  was  matched  to  within  8%.  The  membrane  kinematics  of  the 
finite  element  model  were  similar  to  the  vacuum  test  article  but  overall  membrane 
deflections  predicted  by  the  finite  element  solver  was  less  than  observed  deflections  seen 
in  the  vacuum.  Modal  tuning  brought  the  leading  edge  beam  of  the  wing  to  within  1%  of 
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the  test  article.  Wing  modal  testing  showed  that  the  addition  of  the  membrane  to  the 
beam  reduced  the  first  beam  bending  mode  by  32%  and  the  flapping  motion  was 
dominated  by  the  first  membrane  mode  not  the  first  beam  mode. 
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II.  Literary  Review 


2.1  Low  Reynolds  Number  and  Unsteady  Flow 

As  stated  earlier,  the  small  scale  and  low  Reynolds  number  regime  of  SUASs  coupled 
with  the  unsteady  aerodynamics  associated  with  the  flapping  motion  make  the  flapping 
wing  SUASs  aerodynamics  different  from  conventional  aircraft  control.  Viscous  forces 
become  dominant  for  small  scale  vehicles  and  drag  is  very  sensitive  to  small  changes  in 
lift  [7].  In  conventional  aircraft  aerodynamics,  aerodynamic  coefficients  are  generally 
calculated  or  measured  for  an  aircraft  in  steady  level  flight.  Flapping  wing  are  never  in 
steady  level  flight — they  are  in  a  constant  state  of  motion  in  order  to  maintain  flight  and 
as  such  classical  methods  for  calculating  lift  and  drag  fail  to  add  up  to  measured  values. 
The  quasi-steady  assumption  can  help  simplify  flapping  wing  aerodynamic  calculations 
at  the  expense  of  accuracy.  By  assuming  that  aerodynamic  coefficients  can  be  calculated 
for  the  wing  under  a  series  of  static  conditions,  the  net  aerodynamic  forces  can  be 
calculated  over  the  entire  flapping  cycle  to  get  net  average  values.  This  assumption  is  not 
perfect  and  does  not  take  into  account  unsteady  phenomenon.  As  the  velocity  of  the 
vehicle  increases  the  quasi-steady  contributions  to  net  forces  will  dominate  but 
conversely  as  the  forward  velocity  decreases  the  unsteady  effects  are  greater  [7,  p.  7]. 

The  advance  ratio  is  a  non-dimensional  parameter  of  the  forward  velocity  relative  to  the 
flapping  speed  as  shown  in  Equation  (0.1)  [8,  7,  9]: 


20>  coR 
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where  Vy  is  the  freestream  velocity  of  the  vehicle,  ®  is  the  wing  stroke  amplitude  in 

radians  measured  peak-to-peak,  co  is  the  flapping  frequency  in  Hz,  and  R  is  the  wing 
length.  A  general  rule  of  thumb  states  that  for  J  >  10  quasi-steady  terms  dominate  while 
unsteady  terms  are  the  bigger  players  for  J  <  10  . 

In  2006  Ansari  et  al.  overviewed  the  most  common  ways  in  literature  to  model 
aerodynamics  of  insect  flight  [10].  They  categorized  them  into  the  following  methods: 
steady-state,  quasi-steady,  semi-empirical,  and  fully  unsteady.  First,  they  detailed  the 
main  aerodynamic  phenomena  associated  with  unsteady  effects.  These  aerodynamic 
phenomena  are  the  Wagner  effect,  wake  capture,  apparent  mass  effect,  Kramer  effects, 
and  the  leading  edge  vortex.  These  effects  are  crucial  to  insect  aerodynamics  because  the 
flapping  wings  would  stall  at  the  high  angles  of  attack  seen  in  flight  under  steady-state 
conditions.  This  allows  flapping  wing  vehicles  to  attain  higher  forces  than  those 
predicted  by  quasi-steady  methods. 

The  Wagner  effect  is  an  unsteady  phenomenon  that  actually  reduces  the  effective  lift 
from  quasi-steady  methods.  This  effect,  first  studied  in  1931,  arises  from  the  wing’s 
acceleration  and  deceleration  throughout  the  flapping  cycle  [11].  When  the  wing 
accelerates,  the  circulation  around  it  slowly  rises  to  the  steady-state  equivalent.  During 
this  acceleration,  a  vortex  is  generated  and  shed.  These  vortices  rotate  in  the  opposite 
direction  of  contributory  circulation  and  reduce  the  effective  lift.  However,  Sane 
remarked  that  the  Wagner  effect  may  be  negligible  at  the  Reynolds  number  associated 
with  insect  flight  [1 1].  He  did  acknowledge  that  the  isolation  of  the  Wagner  effect  from 
other  unsteady  phenomena  is  difficult  but  cites  several  recent  studies  that  have  ignored 
the  Wagner  effect. 
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The  second  effect,  wake  capture,  is  also  known  as  wake  passage  and  wing-wake 
interaction.  This  very  significant  effect  was  first  observed  by  Dickinson  in  1994  [11]. 
During  flight,  a  flapping  wing  generates  strong  vortices  and  often  travels  through  its  own 
wake  especially  during  hover.  Strong  vortices  generated  from  stroke  reversal  create 
pockets  of  high  energy  flow.  When  the  wing  reverses  direction  it  interacts  with  these 
high  areas  of  high  energy  flow  and  results  in  a  higher  force  on  the  wing.  Figure  3  below 
pictorially  shows  this  process.  In  2002  Sun  and  Tang  perfonned  a  CFD  analysis  of 
flapping  wing  kinematics  and  concluded  that  the  high  accelerations  associated  with 
stroke  reversal  is  the  cause  of  large  lift  forces  seen  at  the  beginning  and  end  of  a  stroke 
[12].  These  results  contradict  the  experiments  of  Dickinson  in  the  1990s  and  shows  that 
while  the  spikes  in  forces  seen  at  stroke  reversal  are  an  important  part  of  the  flap  cycle, 
different  researchers  attribute  the  cause  for  such  peaks  to  different  phenomenon. 
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Figure  3.  A  hypothesis  for  wing-wake  interaction  from  Sane:  Parts  A-F  depicts  a  wing 
section  as  it  reverses  stroke.  As  the  wing  transitions  from  a  steady  translation  (A)  phase 
and  rotates  around  a  chordwise  axis  in  preparation  for  stroke  reversal,  it  generates 
vorticity  at  both  the  leading  and  trailing  edges  (B).  These  vortices  induce  a  strong 
velocity  field  (dark  blue  arrows)  in  the  intervening  region  (C,D).  As  the  wing  comes  to  a 
halt  and  then  reverses  stroke  (D,E),  it  encounters  this  jet.  As  the  wing  interacts  with  its 
wake,  a  peak  is  registered  in  the  aerodynamic  force  record  (light  blue  arrows),  which  is 
sometimes  called  wake  capture  or  wing-wake  interaction. 
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The  third  effect  is  the  apparent  or  added  mass  effect.  If  an  object  were  to  move  at  a 
constant  velocity  through  an  inviscid  fluid,  kinetic  energy  would  be  conserved  through 
the  separation  of  the  fluid  in  front  of  the  object  to  that  of  it  closing  behind  it.  If  that  same 
object  where  to  accelerate  in  an  inviscid  fluid,  energy  is  not  conserved.  To  account  for 
the  extra  energy  needed  to  accelerate  the  fluid  surrounding  the  object,  a  reaction  force  is 
added  called  the  apparent  or  added  mass  force.  Since  flapping  wings  are  constantly 
accelerating  and  decelerating  this  is  an  important  term.  Sunada  and  Ellington  later 
showed  that  for  wings  with  symmetrical  half-strokes  the  net  apparent  mass  force  is  zero, 
however  not  all  insects  for  flying  vehicles  flap  with  symmetrical  half-strokes  [13]. 

Sane  calls  the  fourth  effect  the  single  greatest  feature  of  insect  wing  aerodynamics: 
leading  edge  vortex.  The  existence  of  leading  edge  vortices  delaying  stall  on 
conventional  aircraft  at  high  angles  of  attack  has  been  studied,  but  direct  evidence  of  their 
role  in  insect  flight  did  not  come  until  Ellington  et  al.  in  1996  [10,  11].  They  built  a  3-D 
scaled  version  of  the  Manduca  sexta  (hawkmoth)  and  seeded  the  flow  field  with  smoke. 

A  leading  edge  vortex  was  observed  through  each  half-stroke  and  was  responsible  for 
delaying  stall  at  high  angles  of  attack.  This  augmentation  on  lift  has  been  described  as 
the  same  phenomenon  seen  on  delta  wings.  Insects  typically  operate  at  high  angles  of 
attack  and  this  effect  is  the  main  contributor.  Since  its  realization,  numerous  studies  have 
been  done  on  the  leading  edge  vortex  for  different  insects.  Many  theories  exist  as  to  why 
the  leading  edge  vortex  is  stable  and  attached  to  the  wing.  Some  studies  have  shown  that 
a  high  spanwise  flow  may  create  this  condition  whereas  others  have  concluded  that  the 
reason  for  its  formulation  changes  with  Reynolds  number  [11],  Although  identified  as  a 
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critical  player  in  flapping  wing  aerodynamics  predicting  leading  edge  vortices  is  very 
difficult  as  they  are  not  fully  understood. 

The  Kramer  effect  is  present  when  a  flapping  wing  rotates  about  a  span-wise  axis 
(change  in  angle  of  attack)  while  translating  through  the  freestream  flow.  When  this 
occurs  the  Kutta  condition  breaks  down  and  additional  circulation  must  be  generated  to 
reestablish  the  Kutta  condition.  Basically,  rotational  circulation  is  generated  to 
counteract  the  change  in  angle  of  attack  of  the  wing.  Depending  on  the  direction  of 
rotation  the  extra  circulation  generated  can  either  supplement  or  detract  from  the  lift  of 
the  wing.  The  Kramer  effect  is  most  prominent  during  stroke  reversal  when  the  wing 
rotates  rapidly. 

Another  aerodynamic  phenomena  detailed  by  Sane  [11]  is  the  clap  and  fling  effect. 
The  clap  and  fling  effect  refers  to  when  an  insect’s  or  any  vehicle’s  wings  touch  or  clap 
at  the  top  of  the  stroke  then  fling  apart.  Theory  behind  the  clap  and  fling  effect  is  that  as 
the  two  wings  clap  together  the  leading  edges  of  the  wings  meet  before  the  trailing  edge. 
As  the  wings  close  together,  the  opposing  circulations  cancel  each  other  out.  The 
vortices  shed  by  the  trailing  edge  are  much  smaller  and  as  a  result  the  Wagner  effect  is 
reduced.  When  the  two  wings  fling  apart  they  create  a  low  pressure  zone  which 
surrounding  air  rushes  to  fill.  This  provides  a  quicker  build  up  of  circulation  or  helps  in 
the  development  of  a  leading  edge  vortex.  The  actual  significance  of  the  clap  and  fling 
effect  is  still  debated.  Many  insects  do  not  use  the  clap  and  fling  mechanism  yet  achieve 
high  lift  values.  Sane  argues  that  in  most  flapping  insects  and  birds  larger  flapping 
amplitude  will  produce  more  lift.  Clap  and  fling,  as  seen  in  nature,  could  simply  be  a 
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result  of  the  insect  trying  to  achieve  the  highest  possible  amplitude  mechanically  possible 
and  there  is  no  added  benefit  to  it. 

These  unsteady  phenomena  are  important  for  a  comprehensive  model  of  a  flapping 
wing.  Full  realization  is  important  to  SUAS  designers  in  order  to  optimize  the  design  for 
specific  mission  capabilities.  The  problem  is  that  these  phenomena  are  extremely 
complicated  and  relatively  recently  discovered  with  most  observations  coming  in  the 
1990s.  As  such,  the  actual  implications  of  these  models  are  not  known  but  researchers 
have  attempted  to  generate  models  to  predict  physical  results.  These  models  can  be  quite 
simple  or  more  complicated  methods  like  CFD  code  that  accounts  for  the  attached  flow  to 
the  wing  as  well  as  the  separated  flow  to  incorporate  unsteady  phenomena  like  the  wake 
capture  effect. 

Steady-state  methods,  as  labeled  by  Ansari  [10],  are  based  upon  momentum  theory. 
The  wings  of  the  insects  are  modeled  as  a  partial  actuator  disc  based  upon  parameters 
such  as  the  vehicle  weight,  wing  length  or  area.  Ansari  notes  that  this  is  a  very  limited 
method  because  it  actually  fails  to  take  into  account  flapping  kinematics  or  wing 
geometry  but  instead  transforms  the  flapping  wing  vehicle  into  a  representative  partial 
actuator  disk. 

Quasi-steady  methods  are  a  fairly  common  way  to  model  aerodynamic  effects  on  a 
flapping  wing.  Osborne  is  generally  cited  as  the  first  contributor  to  quasi-steady  methods 
[14,  9,  10].  One  of  the  most  common  ways  to  model  the  quasi-steady  effects  is  through 
the  use  of  blade  element  theory  (BET).  This  method  is  appealing  because  of  the 
relatively  simple  computational  effort  involved.  Because  of  the  simplifications  of  BET,  it 
generally  under  predicts  lift  forces  because  it  does  not  take  into  account  unsteady 
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phenomena.  This  method  could  be  used  as  a  conservative  estimate  for  SUAS  designers 
early  in  the  design  process.  One  example  of  why  quasi-steady  methods  over  simplify  the 
problem  is  the  small  angle  assumption.  Most  methods  use  thin  airfoil  theory  to  predict 
the  lift  coefficient  on  the  wing.  For  thin  airfoils  at  low  angles  of  attack,  it  was  been 
shown  that  the  lift  coefficient  can  be  approximated  by  CL  =  2na  .  This  assumption  is 

usually  made  for  faster  modeling  purposes  because  it  is  common  for  insect  wings  to 
operate  at  angles  of  attack  of  35°  [10]  which  falls  outside  the  low  angle  of  attack  region. 
However,  quasi-steady  methods  tend  to  be  more  accurate  for  flapping  wing  fliers  at 
higher  speeds.  As  the  speed  increases,  the  unsteady  effects  become  less  significant 
compared  to  the  quasi-steady  terms. 

Semi-empirical  methods  attempt  to  correct  for  the  simplification  of  quasi-steady 
methods  by  supplementing  them  with  corrections  based  on  experimental  results.  These 
methods  can  involve  supplementing  aerodynamic  coefficients  to  try  to  incorporate  the 
unsteady  effects.  Walker  and  Westneat  [15]  attempted  to  added  skin  friction  drag  into 
their  model  based  by  computing  a  skin  friction  coefficient  based  on  Reynolds  number. 
Others  attempt  to  match  lift  and  drag  coefficients  developing  a  unique  equation  for  them 
based  upon  empirically  measured  data.  The  success  of  empirical  methods  is  mixed. 
Empirical  methods  do  build  upon  quasi-steady  methods — increasing  their  accuracy — yet 
the  method  relies  on  experimental  data  rather  than  modeling  the  flow  physics.  As  such, 
the  models  are  very  specific  to  individual  vehicles.  This  method  is  not  entirely  useful  to 
the  aircraft  designer  early  in  the  design  process  as  some  experimental  results  must  be 
attained  to  design  a  vehicle.  Since  these  methods  are  based  upon  quasi-steady 
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assumption,  they  are  not  computationally  intensive  so  if  the  data  is  already  available  then 
the  modeling  can  be  accomplished  relatively  easily. 

Unsteady  methods  are  analytical  methods  that  rely  purely  on  unsteady  aerodynamics. 
The  main  features  of  these  methods  are  that  they  model  the  flow  attached  to  the  wing  as 
well  as  the  wake  from  shedding  leading  and  trailing  edge  vortices.  These  methods  are 
fairly  new  and  as  a  result  a  plethora  of  different  approaches  exist  in  literature.  Some 
methods  can  be  less  computationally  intense  at  the  sacrifice  of  accuracy  while  others 
involve  complex  coding  that  can  be  computationally  intensive. 

BET  was  originally  developed  in  the  late  19th  and  early  20th  century  to  calculate  the 
aerodynamic  forces  on  helicopter  blades  and  is  the  foundation  for  most  helicopter 
aerodynamics  [16,  p.  46].  The  theory  is  based  upon  a  lifting-line  theory  applied  to  a 
helicopter  blade  (modeled  as  an  airfoil)  as  it  rotates  through  the  air.  By  looking  at  the 
cross-section  of  the  blade  as  shown  in  Figure  4,  aerodynamic  coefficients  can  be 
calculated  based  on  the  orientation  and  speed.  For  helicopter  blades  the  total 
aerodynamic  forces  can  be  calculated  by  integrating  the  forces  over  the  entire  blade.  The 
basics  of  this  approach  have  been  adapted  to  predict  flapping  wing  aerodynamic  forces  as 
shown  in  Figure  5. 
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Figure  4.  Example  of  original  blade  element  theory  problem. 


Figure  5.  Blade  element  theory  parameters  definitions. 


Here  </>(t)  is  the  wing  stroke  angle,  R  is  the  wing  length,  c(y)  is  the  chord  length,  and 
dy  is  the  differential  width  of  the  blade  element.  In  this  figure  (/>(l)  is  the  rotation  about 
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the  z-axis  going  into  the  page.  The  lift  on  the  differential  strip  is  derived  from  the  generic 
lift  equation: 


L 


]_ 

2 


pcLv:s 


(0.2) 


dL  =  ^ pCL  (a(t))^2  {t)y2c(y)dy 


(0-3) 


where  L  is  lift,  p  is  air  density,  CL  is  lift  coefficient,  Vx  is  the  freestream  velocity,  and  S 
is  the  wing  planform  area.  In  Equation  (0.3)  it  is  assumed  that  the  only  velocity  is  from 
the  stroke  angular  velocity,  ,  and  dL  represents  the  instantaneous  differential  lift  on 

the  blade  element  section.  The  total  lift  is  found  by  integrating  over  the  wing  so  that  the 
instantaneous  lift  can  be  found  and  similarly  drag  can  be  derived  the  same  way: 


L  =  J dL  =  J \pcl  (« (0)  t1  M  y2c{y)  dy 

0  0  2  (0.4) 

R  R  * 

D  =  \dD  =  \ ~PCD  (a  (t))^2  (t) y2c(y ) dy 

0  0  (0.5) 


where  D  is  the  instantaneous,  and  CD  is  the  coefficient  of  drag.  This  is  a  very  basic 

formulation  of  BET.  In  201 1  Yan  et  al.  developed  and  extended  BET  that  includes  a 
more  comprehensive  model  with  more  parameters  and  unsteady  aerodynamic  methods 
[17]. 

Yan  et  al.  wanted  to  create  a  BET  model  that  captured  the  strong  peak  in  aerodynamic 
forces  seen  at  the  start  of  each  stroke.  They  theorized  that  the  added  mass  effect  and 
wake  capture  were  the  two  unsteady  mechanisms  responsible  for  this  peak;  however, 
since  wake  capture  is  difficult  to  model  they  developed  an  extended  unsteady  blade 
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element  model  that  incorporates  the  added  mass  effect.  They  set  up  their  coordinate 
system  as  shown  in  Figure  8.  Here,  XYZ  is  the  earth  coordinate  system  while  X’Y’Z ’  is 
the  wing-fixed  coordinate  system.  T  represents  the  period  of  a  full  stroke  so  0.5  T  is  one 
half-stroke.  As  before,  (f)(t)  is  the  stroke  angle  and  </>m  is  the  stroke  amplitude.  CLe  and 
CTr  are  defined  as  the  upper  and  lower  limitations  of  the  wing  chord.  Instead  of  using 

along  the  span,  Yan  et  al.  used  a  differential  area  on  the  wing  as  a  blade  element  so  c  is 
the  distance  from  the  axis  of  rotation  (F’),  r  is  the  distance  from  Z’  and  Ac,  Ar  are  the 

differential  widths  of  the  element.  Rb  and  Rr  are  the  upper  and  lower  limitation 
functions  of  spanwise  length  and  are  a  function  of  c  (note  that  Rt  in  Figure  8  should  be 
Rt  ).  With  these  parameters  defined  their  blade  element  model  is  similar  to  the  one  listed 
above  incorporates  more  kinematics  as  seen  below  in  Equation  (0.6)  and  (0.7). 
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Figure  6.  Schematic  from  Yan  et  al.  [17]  detailing  their  formulation  of  extended  blade 

element  theory  model. 
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In  these  equations  F/)an  rot  and  F'ran  rot  are  drag  and  lift  forces  respectively  associated 
with  translation  and  rotation.  The  previous  BET  model  did  not  take  into  account  any 
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translation  of  the  body.  The  model  by  Yan  et  al.  incorporates  this  and  more  through  the 
velocity  term  U  .  The  velocity  term  represents  three  different  parameters: 


(d,c)  +  Utmm  {<j>,r)  +  U bod  str  ($,a) 

(0.8) 

Urot(d,c)  =  cd 

(0.9) 

U,rans  (^)  =  V(j) 

(0.10) 

(y,°0  =  Uhody  cos  (0)  sin  (a) 

(0-11) 

where  Uwt  is  the  instantaneous  linear  component  of  rotational  velocity  of  the  blade 
element,  Utmns  is  the  instantaneous  linear  component  of  stroke  velocity  of  the  element, 

U sir  is  the  instantaneous  flight  velocity  component  in  the  wing  stroke  direction,  and 
Ubody  is  the  flight  velocity  of  the  insect.  In  the  above  equations  a,  a  are  the  angle  of 

attack  and  angular  velocity  of  angle  of  attack  respectively,  (f),  <f)  are  the  wing  stroke  angle, 
and  angular  velocity  of  wing  stroke  angle,  c  is  the  chordwise  distance  from  the  axis  of 
rotation,  and  r  is  the  spanwise  distance  from  the  axis  of  rotation.  The  model  is  more 
comprehensive  because  it  allows  for  further  discretization  of  the  wing  and  incorporates 
body  translational  velocity  as  well  as  the  velocity  associated  with  a  change  in  angle  of 
attack. 

Vanneste  et  al.  [18]  performed  a  similar  approach  on  a  bio-inspired  robotic  insect. 
This  flapping  wing  flyer  was  electro-magnetic  actuated  with  a  wingspan  in  the  3-5  cm 
range.  The  wings  were  micromachined  to  mimic  insect  wings  with  a  high  degree  of 
repeatability.  Vanneste  et  al.  set  out  to  predict  aerodynamic  forces  on  the  wing 
throughout  the  flapping  cycle  through  the  use  of  FE  and  an  aerodynamic  model.  A 
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combination  of  2  node  beam  elements  and  8  node  shell  elements  where  used  in  Ansys.  A 
quasi-steady  approach  was  chosen  for  the  aerodynamic  model  because  the  speed  and  ease 
of  implementation.  Vanneste  et  al.  actually  developed  their  model  after  Sane  et  al.  [19] 
which  is  a  semi-empirical  model.  The  lift  and  drag  coefficients  were  experimentally 
measured  from  a  scaled  fruit  fly  wing  on  a  biomimicking  flapper  immersed  in  an  oil  tank. 
Vanneste  et  al.  successfully  coupled  their  aerodynamic  and  structural  models,  however, 
the  accuracy  of  their  results  cannot  be  determined  as  they  had  no  experimental  results  to 
validate  their  model. 

The  work  of  Seshadri  et  al.  [20]  provided  good  insight  into  the  comparison  of 
classical  steady  methods  versus  unsteady  methods.  They  took  a  rigid  rectangular  wing 
and  hooked  it  up  to  a  four-bar  flapping  mechanism.  To  separate  the  aerodynamic  forces 
from  the  inertial  forces  the  flapping  mechanism  was  placed  in  a  vacuum  and  forces  were 
measured  throughout  the  flap  cycle.  These  forces  were  subtracted  from  the  forces 
measured  by  the  same  flapping  mechanism  flapped  in  air.  Seshadi  et  al.  makes  a  point  to 
mention  that  the  method  of  subtracting  forces  measured  in  a  vacuum  from  the  forces 
measured  in  air  is  only  valid  for  a  rigid  wing.  This  makes  sense  since  a  rigid  wing  will 
roughly  hold  the  same  shape  as  it  flaps  in  a  vacuum  as  it  does  in  air  so  the  instantaneous 
inertial  forces  will  be  the  same  in  both  cases.  For  a  non-rigid  wing,  the  aerodynamic 
forces  will  deform  the  wing  so  that  the  inertial  forces  measured  in  a  vacuum  will  be 
different  from  the  inertial  forces  measured  at  the  same  point  in  the  flap  cycle  in  air. 

After  isolating  the  aerodynamic  forces  the  same  rigid  wing  was  placed  in  a  wind 
tunnel.  Wind  speed  was  varied  to  match  the  range  of  Reynolds  numbers  seen  during  the 
flapping  tests  and  the  angle  of  attack  was  varied  to  find  the  max  possible  lift  coefficient. 
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The  results  from  Seshadri  are  shown  in  Figure  7.  The  wind  tunnel  test  represents  the 
steady  state  method  for  measuring  lift  and  drag  coefficients  based  a  conventional  aircraft. 
The  flapping  wing  tests  contain  all  the  unsteady  effects  not  present  in  conventional 
aircraft.  Unsteady  lift  coefficients  were  up  to  six  times  higher  than  those  measured  by 
steady  state  methods.  This  shows  two  important  points.  One,  flapping  wing  flyers  have  a 
huge  advantage  over  conventional  aircraft  designs  in  the  low  Reynolds  number  regime. 
Two,  the  unsteady  terms  are  very  important  to  aerodynamic  models  to  predict  the  proper 
aerodynamic  coefficients. 


Figure  7.  Experimental  results  from  Seshadri  et  al.  comparing  measured  lift  coefficients 
of  a  rigid  rectangular  wing  in  a  wind  tunnel  (steady)  to  those  measured  on  the  same  wing 
placed  in  a  flapping  mechanism  (unsteady) [20]. 


Seshadri  et  al.  also  developed  an  analytical  model  to  compare  with  the  experimental 
results.  They  used  a  semi-empirical  BET  model  with  aerodynamic  coefficients  based 
upon  a  curve  fit  of  the  unsteady  data  measured  in  the  lab.  While  their  BET  theory  model 
does  use  the  unsteady  aerodynamic  coefficients  measured  in  the  lab  it  does  not  directly 
model  any  unsteady  phenomena  discussed  above  and  instead  tries  to  incorporate  this 
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effects  through  the  use  of  empirical  data.  Their  BET  model  over  predicted  the  lift  forces 
actually  measured  as  seen  in  Figure  8.  Furthermore,  BET  was  only  able  to  capture  two  of 
the  three  peaks  in  lift  over  the  flap  cycle.  Seshadri  et  al.  measured  a  spike  in  lift  during 
stroke  reversal  at  the  end  of  the  first  half-stroke  but  BET  could  not  predict  it.  The  lack  of 
this  third  peak  was  attributed  unsteady  effects  such  as  rotational  circulation  and  wake 
capture.  As  for  the  over  prediction  in  lift  force,  they  hypothesized  that  it  was  due  to  a 
premature  shedding  of  a  leading  edge  vortex  but  had  no  conclusive  data  to  back  up  their 
hypothesis.  Unfortunately  Sehadri  et  al.  did  not  perform  a  BET  model  with  the  steady 
state  coefficients  measured  in  the  wind  tunnel.  It  would  be  expected  that  this  model 
would  under  predict  the  lift,  but  to  what  degree  would  have  been  of  some  interest. 


Figure  37.  Lift  variation  for  (a)  t/’  =  40°  (b)  ifr  =  60°  (c)  if)  =  80° 


Figure  8.  Results  from  Seshadri  et  al.  displaying  their  empirical  BET  model  versus 
experimentally  measured  results  for  stroke  amplitude  of  (a)  1 10°  (b)  100°  (c)  87.3°  [20]. 


In  201 1,  Thielicke  el  al.[21]  created  a  unique  BET  model  incorporating  the  work  on 
delta  wings  into  their  aerodynamic  coefficients.  Leading  edge  vortices  are  a  major  source 
of  lift  in  unsteady  flapping  wing  aerodynamics,  likewise  conventional  delta  wings  also 
rely  on  these  same  vortices  for  lift  especially  at  slow  speeds  during  high  angles  of  attack. 
Thielicke  et  al.  theorizes  that  although  the  stabilization  mechanisms  for  the  leading  edge 
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vortices  may  be  different  from  delta  wings  to  flapping  wings  the  flow  phenomena  and 
aerodynamic  effects  are  similar.  Delta  wings  can  have  lift  coefficients  as  high  as  6  [22] 
which  is  comparable  to  the  increased  lift  coefficients  measured  on  flapping  wings.  They 
use  a  concept  developed  by  Polhamus  of  NASA  back  in  1966  [23]  to  predict  lift 
coefficients  of  sharp  edge  delta  wings.  His  theory,  based  on  potential  flow  and  vortex 
lift,  had  a  trigonometric  relationship  between  aerodynamic  coefficients  and  angle  of 
attack.  Thielicke  et  al.  would  then  use  this  relationship  for  lift  and  drag  coefficients 
within  their  BET  based  model.  The  trigonometric  relationship  from  Polhamus  can  be 
found  in  Equation  (0.12)  below: 

oc 

CL  =K  sin  a  cos2  a +  ^v  cos  a  sin2  a  *^—7  +  CL  (0.12) 

\a  | 

where  CL  is  lift  coefficient,  a  is  angle  of  attack,  CL  is  the  lift  coefficient  at  zero  angle  of 

attack,  and  K  and  Kv  are  the  constants  of  proportionality  in  potential  flow  lift  and  vortex 

lift  respectively.  Drag  coefficient  is  found  by  the  following  relationship  in  Equation 
(0.13): 

CD  =  A CD  +  CDo ,  A CD  =  C,  tan  a  (0.13) 

where  CD  is  the  drag  coefficient,  A CD  is  the  drag  coefficient  due  to  lift,  and  CD  is  the 

drag  coefficient  of  the  wing  at  zero  angle  of  attack.  This  model  would  be  classified  as  a 
semi-empirical  method  since  the  values  of  lift  coefficients  were  based  on  the  research 
done  by  Polhamus  and  the  values  of  CL  and  CD  must  be  experimentally  determined. 

However,  this  is  a  very  innovative  method  compared  to  the  direction  in  recent  literature 
and  the  Thielicke  et  al.  produced  excellent  results  to  back  up  this  approach.  Figure  9 
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shows  an  example  of  the  results  they  obtained  using  their  method.  The  dashed  line 
represents  force  balance  measurements,  the  triangles  represent  the  BET  model  using  the 
aerodynamic  coefficients  from  Polhamus,  and  the  circles  represent  aerodynamic 
coefficients  using  classical  steady  methods.  Their  model  matched  the  experimental 
results  reasonably  well  for  a  semi-empirical  approach.  Figure  10  shows  the  summery  of 
their  results.  The  model  predicted  the  experimental  values  better  than  the  classical 
methods,  but  for  the  horizontal  force  the  method  was  just  as  accurate  and  tended  to  over 
predict  as  the  freestream  velocity  increased.  Thielicke  et  al.  also  noted  that  as  the 
flapping  frequency  increased  the  delta  wing  method  over  predicted  the  forces  generated 
on  the  wing.  This  finding  is  a  cause  of  some  concern  because  the  fastest  frequency  tested 
was  10  Hz.  Insects  such  of  the  fruit  fly  have  a  flapping  frequency  of  200  Hz  [24]  so  if 
the  delta  wing  method  continues  to  over  predict  forces  generated  as  the  frequency 
increases  then  it  would  be  nearly  useless  for  fast  flapping  insect  sized  vehicles. 
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Figure  9.  Results  from  Thielicke  et  al.  showing  the  good  predictability  of  their  “vortex- 
lift”  BET  approach.  The  graphs  show  the  steady  state  coefficients  (circles)  and  “vortex- 
lift”  method  (triangles)  compared  to  force  balance  measurements  (dashed  line).  The 
different  figures  are  for  different  freestream  velocity  (a)  2.28  m/s  (b)  2.57  m/s  (c)  2.84 

m/s.  [21] 
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Figure  10.  Summary  results  from  Thielicke  et  al.  comparing  the  averaged  difference  of 
the  steady  state  coefficient  model  (dark  shade)  and  “vortex-lift”  method  (light  shade) 
from  the  experimental  results  at  different  freestream  velocities.  The  “vortex-lift”  method 
outperformed  the  classical  method  for  the  vertical  force  but  tended  to  over  predict  the 
horizontal  force  as  the  freestream  velocity  increased.  [21] 

The  aerodynamic  forces  associated  with  flapping  wing  flight  are  complex  in  nature 
and  not  yet  fully  understood  since  most  research  in  the  field  began  in  the  1990s.  Many 
models  try  to  capture  these  complex  aerodynamic  interactions.  Some  models  use  simple 
BET  to  achieve  a  fast  approximation  while  others  use  a  combination  of  empirical  terms 
and  unsteady  corrections.  The  highly  unsteady  and  nonlinear  nature  makes  the 
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aerodynamic  forces  difficult  to  predict  but  models  are  getting  more  accurate  as 
researchers  gain  better  understanding  of  the  flapping  wing  effects. 

2.2  Aeroelastic  Response 

Understanding,  predicting,  and  controlling  the  flapping  motion  of  passively 
controlled  wings  is  being  approached  from  many  ways  with  most  of  the  research  focused 
on  the  areas  of  the  aeroelastic  response,  flapping  dynamics,  and  the  structural  qualities 
intrinsic  to  wings  [5,  25].  The  aeroelastic  response  involves  the  interaction  between  the 
wing  structure  and  the  fluid  medium  through  which  it  moves  through.  A  passively 
controlled  wing  deforms  only  by  a  combination  of  the  inertial  forces  of  the  wing  and  the 
aeroelastic  response  so  knowing  and  predicting  the  aeroelastic  response  is  needed  to 
predict  the  behavior  of  a  flapping  wing  small  unmanned  aerial  system. 

Combes  and  Daniel ,  biologists  from  the  University  of  Washington,  performed 
numerous  experiments  with  hawkmoth  wings,  Manduca  Sexta,  in  2003.  They  realized 
that  if  the  aeroelastic  contribution  to  wing  deformation  was  small  compared  to  the 
deformation  caused  by  wing  inertia  alone  then  it  would  be  much  simpler  to  model  insect 
flight  [26].  In  order  to  determine  the  importance  of  inertial  loading  versus  the  fluid 
structure  interaction,  they  sought  to  compare  deformations  of  a  hawkmoth  wing  flapping 
in  air  with  that  in  a  vacuum.  Since  no  fluid  exists  in  a  vacuum,  any  deformation  present 
in  the  wing  will  be  from  inertial  loading  only.  By  comparing  the  deformations  in  a 
vacuum  to  those  in  air,  Combes  and  Daniel  would  be  able  to  detennine  the  dominance  of 
aeroelastic  forces  in  hawkmoth  flight.  Since  they  did  not  have  access  to  a  vacuum 
chamber,  they  created  a  helium  chamber  to  reduce  the  air  density  by  85%  [26],  They  saw 
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slightly  more  pronounced  deformations  in  the  helium  chamber  yet  very  similar  to  the  air 
tested  wings,  leading  to  their  conclusion  that  the  inertial  loading  is  more  significant  in 
predicting  wing  deformation  in  hawkmoths.  To  predict  this  deformation  Combes  and 
Daniel  used  a  simple  finite  element  (FE)  model  of  the  wing  and  adjusted  the  mass 
damping  to  match  the  observed  deformations  in  air  and  helium.  However,  Combes  and 
Daniel  only  picked  three  points  on  the  wing  (wing  tip,  leading  edge,  and  trailing  edge)  to 
match  measure  and  match  to  their  FE  model. 

Norris,  Palazotto,  and  Cobb  set  out  to  corroborate  the  results  of  Combes  and  Daniel 
that  the  hawkmoth’s  flapping  wings  are  ultimately  dominated  by  inertial-elastic  response, 
but  eventually  came  to  the  conclusion  that  the  wing  behaves  more  aeroelastically  [5]. 
Their  approach  differed  from  Combes  and  Daniel  in  two  important  ways:  vacuum  testing 
and  deformation  modeling  [5].  Where  Combes  and  Daniel  were  only  able  to  test  in  a 
helium  chamber,  Norris  et  al  tested  the  hawkmoth  wing  in  a  vacuum  chamber  at  pressure 
of  300  mTorr.  This  test  was  more  representative  of  removing  the  surrounding  fluid  from 
the  wing  and  isolating  the  inertial  deformations  of  the  flapping  wing.  Secondly,  Combes 
and  Daniel  only  modeled  three  points  on  the  hawkmoth  wing  but  Norris  et  al  used 
stroboscopic  photogrammetry  to  record  nearly  the  entire  deformation  throughout  the 
stroke.  Understanding  this  aeroelastic  response  is  important  to  predicting  how  the 
passively  controlled  wing  will  behave. 

As  Combes  and  Daniels  previously  noted  [26],  predicting  the  inertial  forces  is  much 
simpler  than  modeling  the  fluid  structure  interaction  between  the  air  and  the  structure. 
Unlike  conventional  rigid  aircraft  wings,  passively  control  flapping  wings  are  strongly 
coupled  with  the  aerodynamic  forces  produced  throughout  the  flapping  motion.  CFD  is  a 


32 


numerical  method  to  calculate  the  aerodynamic  forces  experienced  by  an  aircraft 
structure  in  flight  commonly  used  nowadays  to  predict  aircraft  response.  Using  CFD  to 
predict  wing  deformations  for  flapping  wing  vehicles  presents  new  and  difficult 
challenges.  When  air  flows  over  the  wing  it  creates  a  pressure  distribution  across  the 
surface  of  the  wing.  For  conventional  rigid  wing  aircraft,  this  pressure  distribution  can  be 
resolved  into  lift  and  drag  forces  that  the  CFD  solver  continues  to  solve  for  consecutive 
time  steps.  For  a  non-rigid  flapping  wing,  the  pressure  distribution  causes  the  wing  to 
deform.  The  deformation  of  the  wing  will  cause  a  change  in  the  resulting  lift  and  drag 
forces.  In  order  to  successfully  model  a  non-rigid  flapping  wing,  the  CFD  solver  must  be 
coupled  with  a  structural  analysis  program.  The  CFD  solver  computes  the  lift  and  drag 
forces  and  sends  these  values  to  a  FE  program  to  calculate  the  wing  deformation  caused 
by  these  forces.  The  defonnation  is  calculated  and  the  new  geometry  is  fed  back  into  the 
CFD  solver.  The  shape  of  the  structure  has  changed  so  the  CFD  solver  must  remesh  or 
adapt  the  aerodynamic  grid  to  fit  this  new  shape.  Remeshing  the  grid  is  a  nontrivial  task 
often  requiring  complex  algorithms  and  multiple  interactions.  This  is  often  the  most 
complicated  and  computationally  intensive  part  of  the  process  as  it  is  done  for  each 
iteration.  The  CFD  software  then  computes  new  lift  and  drag  forces  based  on  this  new 
geometry  and  feeds  back  new  lift  and  drag  forces  and  the  cycle  repeats  as  shown  in  the 
diagram  in  Figure  11.  The  aerodynamic  forces  need  to  be  known  to  calculate  the 
structure  deformation  and  the  structure  shape  needs  to  be  known  to  calculate  the 
aerodynamic  force.  This  process  is  very  computationally  intensive  and  requires  a  fluid 
solver  as  well  as  a  structural  solver. 
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Figure  11.  Schematic  of  workflow  process  between  CFD  and  coupled  FE  solver. 

Nakata  and  Liu  [27]  successfully  coupled  a  FE  based  computational  structural 
dynamic  solver  with  a  CFD  based  dynamic  flight  simulator  to  model  insect  flapping 
flight.  This  fortified  Navier-Stokes  CFD  solver  was  built  in  house  to  specifically  manage 
the  complexities  associated  with  insect  flight  and  the  associated  large  deformation  of  its 
wings.  Nakata  and  Liu  described  how  the  two  meshes  (the  fluid  surface  grid  and 
structure  mesh)  must  deform  together  in  order  to  ensure  no  parts  of  the  load  were  lost 
transferring  from  one  program  to  the  next.  Furthermore,  the  fluid  grid  had  to  adapt 
automatically  as  the  structure  deformed,  requiring  a  sub-iteration  of  the  two  solvers  to 
ensure  proper  convergence  of  the  interaction  between  the  fluid  and  structure  solver  [27]. 
Although  successfully  implemented  by  Nakata  and  Liu,  this  approach  is  very 
computationally  intensive  and  complicated  requiring  specific  solvers  to  handle  the 
specific  case  of  flapping  wing  flight. 
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Lee,  Shin,  and  Lee  used  a  similar  approach  by  using  the  lattice  Boltzmann  method 
(LBM)  and  a  FE  analysis  with  Euler  beam  elements  to  model  structural  deformations  of  a 
flexible  flapping  plate  [28].  The  LBM  is  based  on  the  dynamics  of  particles  used  for 
solving  partial  differential  equations.  The  LBM  essentially  uses  a  fixed  mesh  to  solve 
fluid  structure  interaction  which  is  a  simpler  and  more  efficient  with  a  sacrifice  in 
accuracy  the  greater  the  deformation  [28],  Lee  and  Lee  performed  a  similar  analysis  by 
putting  a  flexible  plate  nonnal  to  free  stream  at  low  Reynolds  numbers  [29].  They  also 
used  the  LBM  and  a  FE  based  solver  to  predict  the  deformations  of  the  thin  flexible  plate 
at  a  range  of  different  Reynolds  numbers. 

A  myriad  of  different  numerical  methods  exists  to  model  a  flapping  wing  with  the 
coupling  of  a  fluid  and  structural  solver  but  they  are  all  computationally  intensive.  The 
goal  of  this  research  is  to  validate  structure  interaction  by  modeling  the  deformations  of  a 
flapping  wing  with  the  aid  of  photogrammetry.  Stereo  photogrammetry  uses  two  images 
taken  of  the  same  surface  to  calculate  a  three  dimensional  location  [30].  With  two 
properly  placed  and  calibrated  cameras,  photogrammetry  can  be  used  to  provide  accurate 
three  dimensional  point  data  of  a  wing  as  it  defonns  throughout  its  stroke.  This  method 
is  especially  appealing  to  for  non-rigid  wings  because  it  is  not  intrusive.  Other  methods 
for  modeling  a  structure  involve  placing  identifying  markers  on  the  structure  in  the  form 
of  a  specific  shape,  paint,  reflective  dots,  etc.  The  application  of  these  materials  to  larger 
rigid  objects  is  negligible  but  the  care  must  be  taken  when  using  those  methods  on 
SUASs.  Because  of  the  light  weight  and  flexible  nature  of  the  wings,  the  weight  and 
stiffness  of  paint  alone  can  be  enough  to  significantly  change  the  flight  mechanics. 
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Several  flapping  wing  studies  have  already  used  photogrammetry  to  achieve  three 
dimensional  models  of  flapping  wings.  Norris  et  al.  used  stroboscopic  photography  and 
the  software  package  PhotoModeler  Scanner  with  only  one  camera  [5],  The  strobe  light 
was  triggered  at  the  same  frequency  as  the  flapping  mechanism  to  “freeze”  the  wing  at 
certain  points  along  the  flapping  stroke.  Image  quality  is  the  upmost  importance  when 
use  photogrammetry  so  they  spent  a  whole  month  perfecting  the  set  up,  lighting,  and 
camera  settings.  They  also  highlighted  that  the  dense  surface  modeling  used  by  the 
software  does  not  work  with  all  objects  but  requires  a  contrasting  or  mottled  surface  [5]. 
Norris  et  al.  used  this  method  on  a  flapping  hawkmoth  wing  both  in  a  vacuum  and  in  air 
to  compare  the  inertial  response  to  the  aeroelastic  response. 

DeLeon  expanded  upon  the  methods  used  by  Norris  et  al  but  instead  of  relying  on 
dense  surface  modeling  decided  to  use  three  dimensional  positions  of  reference  points 
attached  to  the  flapping  wing  [25].  Reference  points  were  used  because  DeLeon  figured 
that  they  allowed  for  more  control  than  dense  surface  modeling  leading  to  more  direct 
comparisons  with  a  higher  confidence  in  the  results.  However  the  presence  of  the 
reference  points  change  the  properties  of  the  wings  so  DeLeon  followed  the  method  of 
Combes  and  Daniel  that  the  application  of  the  points  may  not  exceed  more  than  two 
percent  of  the  weight  of  the  wing  in  order  not  to  significantly  alter  the  natural  frequency 
and  structural  qualities  of  the  wing  [25,  26].  Although  this  process  proved  very  effective 
for  modeling  the  flapping  wing  DeLeon  commented  that  it  was  an  extremely  time 
consuming  process  requiring  manual  inputs  therefore  reducing  the  amount  data  able  to  be 
processed.  This  led  to  the  work  of  Murray  to  improve  upon  the  photogrammetry  data 
capture  process  [30]. 
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Murray  found  that  in  order  to  improve  upon  the  process  that  stereo  photogrammetry 
had  to  be  used  instead  of  the  stroboscopic  photogrammetry.  As  stated  earlier,  stereo 
photogrammetry  uses  two  images  taken  simultaneously  from  cameras  at  different  angles 
to  calculate  three  dimensional  locations  in  the  similar  manner  that  the  brain  uses  the  eyes 
to  compute  depth  perception.  Murray  was  able  to  automate  the  calibration  and  data 
processing  to  allow  for  a  much  less  time  consuming  process.  Murray  confirmed  that 
using  the  Stereo  Vision  Toolbox  produced  reliable  data  and  is  able  to  be  applied  to 
flapping  wing  projects  [30]. 

The  aeroelastic  response  of  a  flapping  wing  involves  the  interaction  between  highly 
complex  aerodynamic  forces  and  the  structure.  While  extensive  research  efforts  have 
gone  into  figuring  out  the  aerodynamic  side  of  the  response,  surprisingly  few  efforts  have 
been  devoted  towards  validating  the  structural  side.  The  inertial  response  of  a  wing  is 
equally  as  important  as  the  aerodynamic  forces  on  the  wing  yet  most  research  never 
attempts  to  validate  their  structural  model.  Most  papers  describe  their  aerodynamic 
process  thoroughly  yet  hardly  mention  the  structural  model  they  use  or  present  any 
results  that  the  structural  model  is  valid.  For  rigid  wings  the  structural  calculations  can 
easily  be  handled  by  a  FE  solver  and  are  hardly  worth  mentioning,  yet  as  research 
explores  the  boundaries  of  flapping  wing  flight,  the  wings  are  becoming  more  and  more 
flexible.  Flexible  wings  are  much  more  difficult  to  model  and  require  a  second  look. 
Flexible  membranes  are  highly  nonlinear  and  one  cannot  assume  that  the  structural  model 
can  correctly  predict  the  response.  The  flexibility  of  wings  can  vary  based  on  the  design 
and  obviously  the  more  rigid  the  wing  the  easier  to  predict  the  response.  This  research 
will  focus  on  a  thin,  highly  flexible  membrane  with  a  modulus  of  elasticity  of 
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approximately  2.5  GPa.  Without  any  stiffening  support  a  membrane  of  this  nature  would 
not  be  able  to  support  forces  normal  to  the  membrane  plane.  This  type  of  membrane  will 
be  very  difficult  to  predict  but  it  is  a  necessary  component  in  predicting  the  aeroelastic 
response  as  more  flexible  wings  are  developed.  This  research  will  set  out  to  validate  the 
deformations  of  a  highly  flexible  wing  model  in  FE. 
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III.  Methodology 


3.1  Test  Article 

The  primary  test  article  was  a  two-wing,  commercial-off-the-shelf  (COTS)  SUAS 
that  can  fly  indoors  and  outdoors  called  the  Duck  Hunter  Extreme  henceforth  referred  to 
as  the  test  article  or  duck.  It  has  a  mass  of  13.5  g,  a  wingspan  of  30  cm  and  a  flap  angle  of 
±15°.  A  picture  of  the  test  article  is  shown  in  Figure  12.  A  remote  control  can  change  the 
flapping  frequency  to  control  thrust  and  the  angle  of  the  horizontal  tale  to  control  the 
flight  direction.  The  test  article  cannot  hover  and  must  maintain  a  forward  velocity  of 
1.75  m/s  in  order  to  achieve  sustained  flight.  Using  Equation  (0.1)  the  advance  ratio  of 
the  wing  can  be  calculated  to  be  0.8  indicating  that  the  test  article  was  well  under  the 
advance  ratio  of  10  and  unsteady  terms  will  dominate  the  quasi-steady  tenns. 

Extensive  measurements  were  taken  of  the  test  article  wing  to  achieve  an  accurate  FE 
model.  The  dimensions  and  thickness  of  two  different  wings  were  taken  in  several 
locations.  The  values  were  averaged  for  use  in  the  FE  model.  The  wing  membrane  was 
composed  of  several  different  types  of  materials.  Most  of  the  membrane  was  composed 
of  a  very  flexible  Mylar,  but  two  strips  of  a  stiffer  material  run  down  the  chordwise 
direction  of  the  wing,  hence  referred  to  as  battens,  as  shown  in  Figure  13.  Along  the 
leading  edge  beam,  a  tape  material  wrapped  around  the  beam  and  attached  to  the 
membrane  on  both  sides.  This  tape  was  the  only  attachment  of  the  membrane  to  the  beam 
and  it  made  the  section  of  the  membrane  adjacent  to  the  leading  edge  stiffer  than  the 
Mylar  but  less  stiff  than  the  battens.  It  is  also  important  to  note  that  there  was  a  rigid 
plastic  tab  located  at  the  rear  of  the  wing  that  connects  the  wing  to  the  SUAS  body  by  a 
pin  joint. 
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The  average  measurements  recorded  are  listed  in  Table  2,  with  their  corresponding 
locations  on  the  wing  shown  in  Figure  13.  The  battens  and  tape  material  were  translucent 
and  difficult  to  see  so  a  corresponding  diagram  of  the  wing  was  constructed  with  the 
stiffeners  and  tape  material  outlined. 


Figure  12.  SUAS  test  article. 
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Table  2.  Wing  Measurements 


Wing  Location 

Part 

Average  Thickness 

(mm) 

1 

Membrane 

0.0215 

2 

Batten  +  Membrane 

0.1426 

3 

Tape  +  Membrane 

0.1185 

4 

Bar 

0.9737 

5 

Bar  +  tape 

1.0613 

6 

Connecter 

0.6172 

7 

Batten  +  Tape  + 

Membrane 

0.2429 

The  duck  achieves  flight  through  a  battery  driven  motor  located  at  the  front  of  the 
fuselage  that  flaps  the  wings.  A  dissected  view  of  the  mechanism  can  be  shown  in  Figure 
14.  This  motor  turned  several  gears  which  in  turn  cranked  a  bar  vertically  up  and  down. 
This  single  bar  was  connected  to  both  the  wing  joints  and  drove  the  flapping  motion  for 
both  wings.  The  joining  of  the  bar  linkage  at  the  top  of  the  duck  had  loose  fittings 
allowing  the  bar  to  slide  during  operation.  Previous  work  on  this  project  had  shown  that 
this  could  be  modeled  as  a  four  bar  linkage  [31].  While  this  mechanism  was  not  a  true 
four  bar  linkage  the  output  to  the  wings  was  relatively  the  same.  This  flapping 
mechanism  was  relatively  simple  when  compared  to  the  complex  motions  found  in 
natural  flyers  like  the  hawkmoth.  The  stroke  plane  angle  and  amplitude  were  fixed  so 
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changing  the  flapping  frequency  was  the  only  available  means  to  modify  the  forces 
generated  during  flight.  This  mechanism  had  its  pros  and  cons.  The  relatively  simple 
nature  did  not  allow  for  much  control  during  flight,  but  replicating  the  flapping 
mechanism  was  simpler.  Since  the  focus  of  this  research  was  on  predicting  the  structural 
response  rather  than  controlling  the  vehicle,  this  made  the  simple  mechanism  an 
appealing  choice. 


Figure  14.  Dissected  view  of  the  mechanism  that  drives  flight  in  the  test  article.  Left 
shows  link  at  the  bottom  of  the  crank  which  corresponds  to  the  top  of  the  upstroke  and 
right  shows  the  bottom  of  the  down  stroke. 


Abaqus  was  chosen  as  the  finite  element  solver  for  this  study  because  it  excelled  at 
handling  highly  nonlinear,  dynamic  problems.  This  strength  was  well-suited  for  this 
research  project  because  most  of  the  wing  was  composed  of  a  flexible  membrane. 

Starting  simple  and  working  toward  a  more  complex  model,  only  the  leading  edge 
beam  was  modeled.  By  starting  simple,  the  model  could  be  easily  validated  at  a  low 
computational  cost,  then  complexity  could  be  added  on.  This  process  was  easy  to  build 
upon  and  debug  when  a  problem  did  occur.  First,  a  beam  only  model  was  created  in 
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NATSTRAN  and  Abaqus.  These  two  models  were  compared  against  each  other  because 
a  relatively  simple  analysis  should  produce  similar  results  whether  from  Abaqus  or 


NASTRAN. 


3.2  Modal  Testing 

A  scanning  laser  vibrometer  was  used  to  perfonn  an  eigenvalue/eigenvector  analysis 
on  the  leading  edge  beam  of  the  duck.  The  beam  of  the  wing  was  clamped  in  between 
two  metal  plates  and  mounted  to  an  electrodynamic  shaker.  The  shaker  would  provide 
the  input,  while  two  lasers  were  used  in  this  experimental  setup  to  measure  the 
input/output  relationship.  A  reference  laser  was  directed  at  the  plate  clamping  the  beam 
to  the  shaker  to  measure  the  velocity  of  the  input  into  the  system.  A  scanning  head  was 
also  set  to  scan  the  length  of  the  beam  to  measure  the  output  velocity.  To  help  strengthen 
the  laser  return  signal,  small  pieces  of  reflective  tape  were  placed  along  the  length  of  the 
beam.  The  tape  was  light  enough  and  the  beam  was  stiff  enough  that  any  effect  on  the 
modal  properties  was  minimal.  A  pseudorandom  input  was  generated  and  input  into  the 
system  through  the  shaker.  Seven  averages  were  recorded  for  each  test  with  800 
frequency  lines  over  a  usable  frequency  range  of  250  Hz  (approximately  80%  of  Nyquist 
frequency).  A  rectangular  window  was  applied  to  the  data  because  the  pseudorandom 
input  was  generated  such  that  it  was  periodic  at  the  Nyquist  frequency.  A  frequency 
response  function  was  generated  from  the  recorded  input/output  signals  using  the  Polytec 
vibrometer  software  and  the  system  eigenvalues  and  eigenvectors  measured  were  used  to 
validate  the  model.  The  test  setup  is  displayed  in  Figure  15  and  a  picture  of  the  beam  in 
the  clamped  configuration  is  shown  in  Figure  15. 
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Figure  15.  Test  setup  for  the  laser  vibrometer.  The  two  lasers  are  located  at  the  bottom 
of  the  image  with  the  scanning  laser  on  the  right  and  the  reference  laser  on  the  left.  In  the 
background  the  wing  is  shown  clamped  to  the  shaker.  The  same  test  setup  was 
performed  for  the  leading  edge  beam. 


Figure  16.  Leading  edge  beam  of  the  duck  wing  shown  clamped  between  two  plates. 
The  plates  were  bolted  to  an  electrodynamic  shaker  for  testing. 


The  FE  model  was  validated  by  comparing  the  first  three  modes  to  the  measured 
modes  in  the  laboratory.  Any  discrepancies  between  the  model  and  the  experimental 
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results  were  reconciled  via  model  tuning.  The  density  of  the  model  was  verified  by 
weighing  the  beam  and  changing  the  FE  model  density  to  match  the  measured  weight.  A 
clamped  beam  is  easy  to  model  so  a  representative  model  was  constructed  in  MATLAB. 
An  eigenvalue  analysis  was  run  to  ensure  the  modes  matched  the  Abaqus  model. 
Modeling  tuning  occurred  by  changing  the  modulus  of  elasticity  until  the  first  two  modes 
matched  the  experimental  results.  A  sequential  quadratic  programming  optimization 
algorithm  in  MATLAB  modified  the  modulus  of  elasticity  until  the  first  two  modes 
matched.  The  model  had  been  tuned  to  the  clamped  beam,  but  the  housing  within  the 
duck  that  the  leading  edge  beam  sits  in  has  some  wobble.  This  boundary  condition  was  a 
source  of  error  between  the  model  and  the  physical  test  article. 

To  remedy  this  unknown  factor,  modal  testing  was  accomplished  with  the  entire  test 
article.  The  duck  was  strapped  to  a  plate  via  zip  ties  and  the  plate  was  bolted  to  the 
shaker.  The  Styrofoam  body  of  the  duck  did  not  allow  direct  mounting  to  the  shaker  but 
the  zip  ties  firmly  secured  the  duck  so  not  to  allow  the  body  to  slip  or  wiggle  during 
testing.  The  bar  was  mounted  in  housing  on  the  duck  in  place  of  the  wing  (Figure  17) 
and  the  scanning  laser  vibrometer  was  used  in  the  same  manner  as  before.  This  gave  a 
more  realistic  boundary  condition  and  allowed  for  insight  into  how  the  modes  were 
affected  by  the  compliance  in  the  joint. 
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•  •  •  • 

•  •  •  • 


Figure  17.  Top  view  of  the  test  article  used  for  modal  testing  with  the  wing  membrane 
removed  from  the  leading  edge  beam  on  the  left  side.  This  configuration  would  show 
different  results  than  the  bar  in  the  clamped  configuration. 

The  model  was  then  tuned  by  modifying  the  boundary  condition.  The  clamped 
boundary  condition  was  stiffer  than  the  actual  joint.  To  weaken  the  attachment  point  the 
modulus  of  elasticity  and  length  of  the  element  attached  to  the  boundary  condition  was 
varied  until  the  model  matched  the  first  two  measured  modes.  The  methodology  behind 
this  decision  was  that  the  full  bar  and  wing  model  would  be  actuated  by  forcing  a  rotation 
at  the  root  of  the  wing.  In  the  Abaqus  solver  this  can  only  be  accomplished  by 
constraining  the  degrees  of  freedom  at  that  point  and  then  applying  the  forced  rotation. 

By  weakening  the  beam  at  the  attachment  point,  it  would  simulate  this  weak  connection 
and  the  results  could  be  properly  seen  in  the  dynamic  analysis  as  well  as  the  frequency 
analysis. 

Again  MATLAB  was  used  to  optimize  this  connection.  This  optimization  gave  an 
accurate  bar  model  with  an  appropriate  boundary  condition.  With  this  simple  model 
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validated,  the  next  step  was  to  add  the  wing  membrane  to  the  leading  edge  bar.  As  with 
the  bar,  a  wing  model  was  created  in  Abaqus  and  FEMAP.  Analyzing  these  models  gave 
a  comparison  between  the  NASTRAN  and  Abaqus  solver. 

Two  wings  were  measured  in  detail  using  a  micrometer  for  the  thickness  of  the  wing 
membrane,  batten,  connector,  tape,  tape/membrane  overlap,  and  the  bar  thickness.  Figure 
13  above  shows  the  locations  for  all  the  measurements  and  Table  1  lists  the  average 
measurement  of  each  part.  Next,  material  property  testing  was  conducted  since  the  actual 
materials  used  for  each  part  was  unknown.  Several  spare  wings  were  dissected  into  their 
different  components:  leading  edge  beam,  wing  membrane,  batten,  and  tape.  The 
connector  was  not  tested  because  it  was  fairly  stiff  and  its  affect  on  the  outcome  of  the 
simulation  was  insignificant.  The  tape  and  membrane  overlap  area  was  too  small  to  test 
in  the  appropriate  testing  equipment. 

To  model  the  membrane  portion  of  the  wing  more  accurately  in  Abaqus,  a  picture  of 
the  wing  was  taken  from  directly  above.  This  image  was  loaded  into  Abaqus  and  the  part 
was  crafted  by  tracing  the  wing  outline  and  all  other  features.  The  leading  edge  was 
modeled  as  a  simple  beam  as  a  separate  part  from  the  rest  of  the  wing.  The  membrane 
was  sectioned  so  that  each  feature  of  the  wing  could  be  properly  segregated.  The  final 
assembly  can  be  found  in  Figure  18.  More  details  about  the  model  can  be  found  below. 
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Figure  18.  The  Abaqus  wing  model  includes  several  features  found  on  the  duck  such  as 
the  tape,  battens,  and  connector.  Note:  the  leading  edge  beam  runs  the  whole  length  of 
the  wing  membrane.  The  wing  membrane  is  connected  to  the  beam  by  the  strip  of  tape  at 

the  leading  edge. 


With  the  full  model  developed,  modal  analysis  was  then  performed  on  the  entire  wing 
with  methodology  similar  to  that  used  for  the  bar.  The  setup  for  the  wing  testing  was  the 
same  as  shown  in  Figure  15  above.  The  wing  had  to  be  coated  with  Spotcheck  for  the 
scanning  vibrometer  to  receive  a  strong  enough  return  from  the  laser.  Spotcheck  is  a 
spray  of  white  developing  particles  that  can  be  quickly  applied.  Without  this  layer  of 
particles,  most  of  the  scanning  laser  would  pass  through  the  translucent  membrane.  This 
method  added  some  weight  to  the  wing,  which  was  taken  into  account  when  model 
verification  occurred,  but  was  a  necessary  step  to  accomplish  this  testing.  The  extra 
weight  of  the  Spotcheck  most  likely  made  the  measured  modes  lower  than  those  seen 
during  flight.  Spotcheck  was  applied  sparingly  to  reduce  the  effect  of  the  added  weight 
but  it  was  necessary  to  achieve  a  sufficient  return  from  the  laser.  The  wing  was  tested  in 
a  clamped  configuration  (Figure  19)  as  well  as  in  the  configuration  used  in  flight  (Figure 
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20).  The  modal  modeling  was  accomplished  in  Abaqus.  As  previously  stated  with  the 
beam,  the  test  article  housing  had  some  compliance  and  was  not  a  perfect  clamped 
boundary  condition  so  the  modal  tuning  was  accomplished  in  the  same  manner  as 
described  for  the  beam. 


Figure  19.  Duck  wing  attached  to  shaker  in  clamped  configuration.  The  wing  was 
painted  with  Spotcheck  to  achieve  a  strong  laser  return.  The  reference  laser  can  be  seen 
on  the  clamp  and  the  scanning  laser  on  a  point  in  the  middle  of  the  wing. 


50 


Figure  20.  Duck  clamped  to  shaker  for  modal  testing. 


3.3  Air  and  Vacuum  Testing 

In  conjunction  with  modeling  the  wing,  the  test  article  was  also  tested  in  air  and 
vacuum.  Predicting  the  aeroelastic  response  of  the  wing  is  difficult  because  of  the 
complex  interaction  between  the  aerodynamic  and  inertial  forces.  As  discussed  earlier 
the  aerodynamic  forces  themselves  are  difficult  to  understand,  let  alone  model.  When 
they  are  coupled  with  the  inertial  response  the  two  forces  become  difficult  to  separate. 

For  a  rigid  wing  the  inertial  forces  are  easier  to  isolate  because  aerodynamic  forces  have 
little  effect  on  the  wing  shape  as  it  flaps.  Hence,  the  inertial  forces  can  be  found  by 
testing  the  wing  in  a  vacuum  since,  in  the  absence  of  air,  the  only  forces  measured  will  be 
inertial  forces  as  detailed  by  [20]. 

The  flexible  wing  presents  more  difficulties  because  the  wing  will  deform  in  a 
different  manner  in  a  vacuum  than  under  the  influence  of  aerodynamic  forces.  To  test 
this  assumption  the  wing  was  tested  in  a  vacuum  as  well  as  in  air.  The  vacuum  testing 
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was  important  because  it  allowed  for  a  qualitative  assessment  of  the  membrane  deflection 
in  air.  As  discussed  earlier  the  FE  model  does  not  incorporate  aerodynamics  so  the 
deflections  seen  in  the  model  should  compare  well  with  deflections  measured  in  a 
vacuum. 

Vacuum  testing  occurred  in  a  custom  vacuum  chamber  at  AFIT  that  allowed  the 
pressure  to  range  from  atmospheric  to  less  than  one  Torr.  This  vacuum  chamber  was 
ideal  because  all  the  sides  were  translucent,  which  allowed  for  easy  viewing  of  the  duck 
from  all  sides.  Three  high-speed  cameras  were  setup  to  capture  the  wing  throughout  the 
flap  cycle.  This  setup  would  allow  for  qualitative  assessment  of  the  membrane  billow 
and  a  quantitative  assessment  through  the  use  of  photogrammetry.  A  picture  of  the  setup 
can  be  found  in  Figure  2 1 .  Two  of  the  cameras  were  setup  on  either  side  of  the  vacuum 
chamber  so  that  one  looked  at  the  side  of  the  duck  while  the  other  captured  the  front  of 
the  duck.  The  third  camera  was  placed  directly  above  the  duck.  Synchronization 
between  all  three  cameras  was  accomplished  through  a  trigger  controlled  by  a  MATLAB 
program. 
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Figure  2 1 .  Experimental  setup  of  vacuum  testing.  The  test  article  was  placed  in  the 
vacuum  chamber  and  wing  stroke  was  capture  by  three  high  speed  cameras  (right,  bottom 

left,  and  top  of  photo). 


Under  normal  operating  conditions  the  duck  was  powered  by  a  rechargeable  lithium- 
ion  battery.  Placing  this  battery  in  a  vacuum  could  be  potentially  hazardous.  To 
circumvent  this  possibility  the  battery  was  removed  and  the  duck  hooked  up  to  an 
equivalent  DC  voltage  power  source  outside  the  vacuum  chamber.  With  this  setup,  the 
duck  was  safely  controlled  outside  the  vacuum  chamber  with  the  use  of  the  remote 
control.  The  high-speed  images  would  give  a  qualitative  analysis  of  the  membrane  while 
photogrammetry  could  provide  a  3D  model  of  the  wing  deflection  as  a  function  of  time. 
Air  testing  occurred  in  two  different  manners.  One  occurred  while  the  duck  was  in  the 
vacuum  chamber  under  atmospheric  conditions — hooked  up  to  the  DC  power  supply.  A 
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second  test  occurred  with  another  test  article  flapping  with  the  commercially  supplied 
battery.  The  difference  in  these  two  tests  would  reveal  any  difference  in  the  two  different 
power  supplies. 

3.4  Structural  Model 

As  stated  earlier  the  primary  development  of  the  structural  model  was  performed  in 
Abaqus.  The  final  model  was  composed  of  two  parts:  the  leading  edge  beam  and  the 
wing  membrane.  The  final  material  properties  detennined  through  testing  can  be  found 
in  Table  3.  These  material  properties  represent  the  base  model  in  Abaqus.  A  multitude 
of  different  models  were  run  changing  various  parameters  to  see  their  effect  on  the 
output.  This  was  an  effort  to  match  the  physical  deformations  seen.  Membrane  elements 
and  plate  elements  were  both  used  in  the  various  models  and  the  beam  was  always 
modeled  with  beam  elements.  Since  membrane  elements  only  resist  in  place  forces,  they 
offered  a  lower  computational  cost  with  the  risk  of  poor  accuracy.  The  plate  elements 
provide  more  degrees  of  freedom  at  a  longer  compute  time.  Various  combinations  of 
these  were  tried  to  optimize  the  ratio  between  accuracy  and  computational  intensity. 
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Table  3.  FE  Model  Properties 


FE  Spec 

Beam 

Membrane 

Batten 

Connector 

Modulus  of 

Elasticity  (GPa) 

85.9 

2.47 

5 

50 

Poisson’s  Ratio 

0.05 

0.489 

0.489 

0.05 

Density  (kg/m  ) 

1584 

917 

917 

917 

Element  Type 

Beam 

Membrane 

Shell 

Shell 

The  leading  edge  beam  and  the  wing  membrane  were  connected  in  Abaqus  through  a 
tie  constraint  along  the  nodes  of  the  beam  allowing  the  assembly  to  act  as  one  part.  The 
only  load  acting  on  the  structure  was  gravity.  A  pinned  boundary  condition  was  created 
at  the  inner  circle  on  the  connector  to  simulate  the  connection  to  the  test  article  body. 
Since  the  wing  on  the  test  article  was  actuated  by  a  four  bar  mechanism  described  above 
a  forced  rotation  was  placed  along  the  part  of  the  leading  edge  beam  jutting  out  from  the 
wing  membrane.  Many  other  methods  were  attempted  to  achieve  the  proper  flapping 
motion  (force  at  the  root,  moment  input,  etc.)  but  a  forced  rotation  was  the  only  method 
that  produced  repeatable  and  reliable  results.  Because  of  the  complexities  associated 
with  dynamic  simulations  two  different  meshes  were  created.  One  was  a  coarse  mesh 
used  for  most  of  the  analysis  because  it  allowed  for  a  reasonable  visualization  of  the 
results  at  a  lower  compute  time.  If  more  exact  data  was  needed  or  once  a  model  was 
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verified  a  finer  mesh  would  be  used.  Figure  22  and  Figure  23  below  show  an  example  of 
a  coarse  mesh  and  fine  mesh  used  in  the  analysis. 


Figure  22.  An  example  of  a  coarse  mesh  used  for  most  analysis. 


Figure  23.  An  example  of  a  finer  meshed  used  for  in-depth  analysis. 

The  dynamic  nature  of  the  flapping  motion  made  the  structural  problem  difficult  to 
solve.  Abaqus  has  two  different  solvers  for  dynamic  simulations:  implicit  and  explicit. 
In  an  implicit  dynamic  analysis  the  entire  set  of  nonlinear  equilibrium  equations  must  be 
solved  at  each  time  increment.  The  implicit  analysis  could  take  a  long  time  but  the 
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equations  are  solved  at  each  time  step  making  this  method  unconditionally  stable.  For 
example  there  is  no  limit  on  the  size  of  the  time  increment  for  an  implicit  analysis.  The 
explicit  analysis  does  not  solve  the  nonlinear  equations  but  rather  used  known  quantities 
at  previous  time  steps  to  estimate  the  displacements  and  velocities  at  a  new  time 
increment.  This  method  was  generally  less  computationally  intensive  but  the  time  step  is 
conditionally  stable.  There  exists  a  limit  on  the  size  of  the  time  step  before  the  solution 
becomes  unstable.  Abaqus  defines  this  limit  as  the  time  it  takes  a  stress  wave  to  cross  the 
smallest  element  length  in  the  model.  The  smaller  the  element  the  shorter  the  time  step 
can  be  taken.  Another  downside  of  the  Abaqus  explicit  solver  was  that  fewer  element 
types  are  supported  for  this  analysis.  However,  beam,  plate,  and  membrane  elements 
were  all  supported  so  that  both  the  implicit  and  explicit  solvers  could  be  tested  against 
each  other. 
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IV.  Results 


4.1  Modeling  Goal 

The  purpose  of  this  research  was  to  correctly  model  the  structural  response  of  a 
flapping  wing  of  the  proposed  test  article.  To  that  end,  before  modeling  was  started  it 
was  beneficial  to  test  the  actual  model  in  air  so  that  the  wing  kinematics  could  be 
observed.  Before  quantitative  data  was  taken,  high  speed  images  of  the  wing  gave  a 
frame  of  reference  for  the  modeling  that  would  occur  later.  The  results  from  the  model 
were  qualitatively  compared  against  the  images  to  provide  a  gauge  to  judge  if  the  model 
was  generally  predicting  the  motion  or  if  something  was  not  correct.  This  comparison 
saved  time  up  front  because  it  allowed  for  a  quick  comparison  to  the  test  article.  Even 
though  the  FE  results  did  not  model  air  and  therefore  should  not  match  the  wing  in  air,  it 
was  expected  that  the  kinematics  of  the  wing  to  be  similar.  It  was  theorized  that  while 
the  vacuum  results  would  probably  show  greater  deflections,  the  overall  movement  of  the 
wing  would  be  the  same. 

The  high  speed  cameras  were  used  to  capture  images  of  the  flapping  wing  as  it  was 
clamped  to  a  test  stand.  The  series  of  images  in  Figure  24  show  the  motions  of  the  wing 
throughout  the  flap  cycle.  At  the  beginning  of  the  down  stroke,  the  wing  beam  began  to 
move  down  while  most  of  the  membrane  movement  was  still  in  the  direction  of  the 
preceding  stroke.  The  wing  membrane  underwent  large  deformations  during  the  stroke 
with  the  most  prominent  during  stroke  reversal.  During  the  down  stroke  the  beam  again 
led  membrane  until  the  beam  was  about  level  then  the  aft  part  of  the  membrane  surpassed 
beam.  When  the  beam  began  the  upstroke,  most  of  the  membrane  was  still  moving  down 
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resulting  in  a  larger  defonnation.  At  the  end  of  the  upstroke,  the  beam  arrived  at  the  top 
just  before  the  aft  part  of  the  membrane. 
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(d)  75%  of  Stroke 


(e)  3%  of  Stroke 

Figure  24.  The  flap  cycle  synchronized  to  percent  of  the  stroke  period  [31] 
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Figure  25  shows  a  trace  of  the  tip  displacement  during  multiple  strokes.  Birds  and 
insects  with  complex  flapping  kinematics  can  have  unique  wing  tip  traces  as  shown  in 
Figure  26[24].  The  test  articles  wingtip  path  was  relatively  simple  but  interestingly 
enough,  the  wingtip  path  most  closely  resembled  a  fruit  fly  (two  wings)  from  the  diagram 
below  whose  cruising  speed  is  3  m/s  (the  test  article’s  cruising  speed  is  1.75  m/s).  The 
wingtip  trace  for  the  test  article’s  wing  was  consistent  over  multiple  flap  cycles  and 
simple,  which  was  another  reason  this  test  article  makes  a  good  choice  for  FE  modeling. 
The  features  observed  during  this  testing,  specifically  the  leading  of  the  beam  or 
membrane  during  different  phases  of  the  stroke,  wingtip  path,  and  consistent  deflections 
during  each  flap  cycle,  were  used  as  qualitative  gauges  of  accuracy  throughout  the 
modeling  process. 
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(a)  Side  View  (b)  Front  View 

Figure  25.  The  path  of  the  tip  can  be  consistently  measured  from  the  images  [31]. 
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Figure  26.  Wingtip  trace  for  (a)  albatross  fast  gait,  (b)  pigeon  slow  gait,  (c)  horseshoe  bat 
fast  gait,  (d)  horseshoe  bat  slow  gait,  (e)  blowfly,  (f)  locust,  (g)  June  beetle,  (h)  fruit  fly 
[24].  The  test  article  most  closely  resembles  the  fruit  fly  wingtip  trace. 


4.2  FE  Solver 

Starting  with  the  least  complex  system,  the  dynamics  of  the  leading  edge  beam  were 
modeled  in  NASTRAN  and  Abaqus.  While  NASTRAN  is  a  proven  finite  element  solver 
it  was  questionable  how  well  it  would  handle  the  dynamic  and  highly  nonlinear 
interactions  of  a  wing  membrane.  Abaqus  on  the  other  hand,  is  known  for  handling 
dynamic  problems,  but  was  a  relatively  newer  program.  Modeling  only  the  beam  was  a 
first  step  towards  modeling  the  entire  wing  but  it  was  also  gave  early  insight  as  to  which 
solver  would  better  handle  the  nonlinear  dynamic  analysis.  These  early  models  were 
developed  with  the  material  properties  used  in  the  previous  work  on  the  duck  [31] 
because  material  testing  had  not  yet  been  accomplished.  One  parameter  that  will  be  used 
constantly  throughout  the  course  of  this  paper  to  compare  different  models  is  the  tip 
deflection  over  time.  The  tip  deflection  is  simply  a  trace  of  the  tip  of  the  end  of  the 
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leading  edge  beam  during  the  flapping  cycle  as  shown  in  Figure  25.  It  was  theorized  that 
if  a  model  could  reproduce  the  tip  displacement  over  time  for  the  test  article,  then  the  rest 
of  the  beam  would  match  closely.  Tip  deflection  gave  a  quantifiable  parameter  that  could 
easily  be  obtained  in  FE  models  and  was  measurable  in  the  laboratory. 

Figure  27  below  shows  the  tip  deflection  versus  time  for  the  various  models.  All 
models  were  started  at  the  bottom  of  the  stroke  because  the  velocity  at  this  point  would 
be  zero  and  this  would  allow  for  a  quicker  damping  out  of  transients  than  if  started  during 
midstroke.  For  a  frame  of  reference,  the  rigid  body  motion  of  the  beam  has  been  plotting 
on  the  graph.  This  series  would  be  the  tip  deflection  of  the  leading  edge  beam  if  it  were 
to  stay  perfectly  rigid  without  any  deflection  throughout  the  flapping  cycle.  The 
NASTRAN  methods  both  reached  a  steady  state  around  0.3  seconds  and  both  showed 
very  similar  deformations.  However,  the  steady  state  achieved  was  not  purely  sinusoidal. 
Higher  order  modes  seemed  to  be  present  so  that,  after  steady  state  was  achieved,  every 
other  down  stroke  was  larger  than  the  average  amplitude. 

This  analysis  did  not  shed  any  light  into  the  accuracy  of  one  method  over  the  next 
because  of  the  absence  of  experimental  results  for  comparison.  Both  methods  predicted 
approximately  a  1 10  mm  max  deflection  (max  to  min  peak)  at  some  point  during  the  flap 
cycle  while  the  rigid  body  calculations  only  showed  a  deflection  amplitude  of  70  mm. 
Roughly  30  mm  of  deflection  was  caused  by  bending  of  the  beam  at  stroke  reversal  but 
there  was  also  the  existence  of  two  peaks  at  the  height  of  the  upstroke.  The  Abaqus 
explicit  and  implicit  solvers  both  predict  almost  identical  results  at  the  beginning  of  the 
stroke;  yet  over  time  the  implicit  solver  grew  with  each  period  until  it  seemed  that  the 
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solution  was  growing  unbounded.  This  finding  suggested  that  the  implicit  analysis  was 
unstable  so  this  method  was  investigated  further. 


Figure  27.  Tip  displacement  over  time  comparison  between  NASTRAN  and  Abaqus  of 

the  leading  edge  beam  only. 

While  the  implicit  solver  seemed  to  perform  well  0.6  seconds  into  the  analysis,  after 
this  point  the  amplitude  began  to  grow  with  each  stroke  until  it  seemed  the  solution  was 
growing  unstable.  This  trend  was  disconcerting  because  although  the  implicit  solver  was 
supposed  to  be  more  computationally  intensive,  it  was  supposed  to  be  the  more  stable  of 
the  two  solvers.  The  same  model  was  run  reducing  the  step  size  from  0.001  to  0.0001  to 
see  the  effect  on  the  solution.  Changing  the  step  size  had  a  drastic  effect  on  the  solution 
as  seen  in  Figure  28.  While  the  solution  was  relatively  the  same  through  0.4  seconds,  the 
implicit  analysis  with  a  larger  time  step  increased  over  time  especially  on  the  bottom  of 


64 


the  down  stroke.  This  gave  less  confidence  in  the  implicit  solver,  because  it  seemed  that 
the  solver  could  not  obtain  convergence  with  that  large  of  a  time  step.  Coupled  with  the 
fact  that  the  implicit  solver  took  much  longer  to  run  at  a  time  step  of  0.0001,  the  explicit 
analysis  was  chosen  as  the  solver  of  choice  because  the  addition  of  the  wing  membrane 
would  only  increase  the  complexity  and  analysis  time  drastically,  and  for  this  scenario  the 
explicit  solver  seemed  to  handle  the  nonlinearities  better  than  the  implicit  solver  did. 
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Tip  Displacement  (mm)  Tip  Displacement  (mm) 


Figure  28.  Comparison  of  implicit  solvers. 


Investigation  of  the  best  solver  to  use  began  with  incorporating  the  full  wing  into  the 
NASTRAN  solver.  The  same  geometry  in  the  model  developed  by  Alerding  [31]  was 
imported  into  FEMAP.  The  beam  and  wing  model  in  FEMAP  is  shown  in  Figure  29. 
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The  model  was  adequately  meshed  and  analysis  was  perfonned  using  the  NASTRAN 
linear  transient  and  nonlinear  solvers.  Beam  elements  were  used  to  mesh  the  beam  and 
plate  elements  that  comprised  the  membrane.  Figure  30  shows  the  membrane  deflection 
predicted  by  these  analyses.  The  membrane  deflections  were  completely  unreasonable, 
characterized  by  a  high  number  of  ripples  deflecting  unrealistically  large  amplitudes. 

The  membranes  looked  as  if  they  exhibited  an  extremely  low  stiffness.  Several  attempts 
were  made  to  remedy  these  results  in  the  model  but  the  only  way  to  force  the  results  to 
qualitatively  match  reality  was  to  arbitrarily  add  bending  stiffness  to  the  membrane.  It 
was  determined  that  the  NASTRAN  solver  could  not  handle  the  highly  nonlinear  aspects 
of  the  wing  membrane  so  Abaqus  was  chosen  as  the  FE  solver. 


Figure  29.  Wing  modeled  and  meshed  in  FEMAP  for  comparison  of  the  NASTRAN 

with  the  Abaqus  solver. 


67 


Figure  30.  Membrane  mesh  deformations  are  unrealistic  and  fail  to  match  physical 
model.  Left:  NASTRAN  linear  transient  solver.  Right:  NASTRAN  advance  nonlinear 

solver. 


4.3  Material  Testing 

Material  testing  on  the  wing  produced  some  mixed  results.  Standard  tensile  testing 
was  performed  on  the  membrane,  battens,  and  tape.  The  wing  membrane  gave  repeatable 
results  with  the  average  modulus  of  elasticity  at  2.47  GPa.  This  find  also  corroborates 
well  with  the  work  of  Gogulapati  [32]  who  tested  a  similar  wing  membrane  material  and 
found  a  modulus  of  elasticity  was  2.72  GPa.  The  tape  results  were  unreliable  due  to  the 
adhesive  between  the  layers.  The  tape  made  this  material  behave  as  a  composite  with  a 
very  weak  adhesive.  During  testing  the  material  would  stretch  and  then  part  of  the 
adhesive  would  give,  then  the  material  would  stretch  more  before  more  adhesive  would 
give.  The  adhesive  skewed  the  stress  strain  curves  so  in  the  linear  region  it  appeared  that 
the  tape  and  membrane  layer  had  an  average  modulus  of  elasticity  of  0.839.  This  was 
deemed  unrealistically  low.  For  the  final  model  the  membrane  modulus  of  elasticity  was 
used  for  the  tape  area. 
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The  batten  material  was  much  stiffer  than  the  membrane.  As  a  result  the  load  cell  on 


the  testing  apparatus  was  maxed  out  (40  N)  before  the  material  failed.  Unfortunately  the 
next  load  cell  size  up  would  not  have  been  able  to  accurately  measure  the  small  loads 
associated  with  this  material.  For  the  model  a  conservative  estimate  of  5  GPa — 
approximately  twice  the  membrane  modulus — was  used.  For  reasons  stated  above  the 
load  cell  could  not  measure  the  relatively  stiff  carbon  fiber  beam.  As  a  work  around  the 
beam  was  divided  into  several  sections  and  tested  in  a  three  point  bending  configuration 
shown  below  in  Figure  3 1 .  This  three  point  bending  gave  repeatable  results  and  a 
modulus  of  elasticity  of  88. 1  GPa  was  found  for  the  leading  edge  beam.  The  same  test 
setup  used  for  the  leading  edge  beam  was  also  used  to  test  the  battens  in  an  attempt  to 
find  the  material  properties.  Unfortunately,  the  battens  did  not  exhibit  enough  bending 
stiffness  so  the  forces  generated  were  below  the  noise  threshold  for  the  sensors.  As  part 
of  the  material  testing,  the  leading  edge  beam  was  weighed  and  found  to  have  a  density 
of  1584  kg/m  .  The  density  of  the  material  was  updated  in  the  Abaqus  model  to  reflect 
these  material  property  changes. 


Figure  3 1 .  Three  point  testing  configuration  for  material  testing.  On  the  right  the 
apparatus  is  in  the  open  configuration.  On  the  left  the  material  (batten  in  this  picture)  is 

bent. 
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4.4  Laser  Vibrometer 


Model  tuning  via  modal  testing,  as  outlined  in  Chapter  3,  proved  very  insightful.  The 
frequency  response  functions  (FRFs)  and  coherence  for  the  leading  edge  beam  in  the 
clamped  configuration  are  shown  below  in  Figure  32.  The  coherence  was  excellent  over 
the  range  of  interested  frequencies.  The  first  three  modes  are  shown  in  Table  4.  The 
Abaqus  bar  only  model  in  a  perfectly  clamped  configuration  gave  similar  results  but  there 
was  still  a  1-3%  error  in  each  mode.  To  tune  the  model,  a  MATLAB  FE  code  was 
developed  and  ran  to  show  that  they  matched  the  Abaqus  results.  As  discussed  in 
Chapter  3  an  interior-point  algorithm  was  used  for  all  optimization  routines.  The 
modulus  of  elasticity  was  modified  to  match  the  first  mode  in  the  beam.  The  results  of 
the  optimization  are  shown  in  Table  4.  The  model  was  tuned  within  1%  of  the  physical 
beam  and  the  modulus  of  elasticity  was  decreased  by  2.5%  to  achieve  this.  The  model 
was  tuned  to  match  the  first  mode  because  this  mode  occurred  nearest  the  flight  flapping 
frequency  of  15  Hz.  Note  that  while  the  first  mode  matches  by  design  the  second  and 
third  modes  still  exhibit  some  error.  The  two  plates  the  bar  was  sandwiched  between  was 
not  a  perfect  clamp  so  some  error  was  expected. 

The  leading  edge  beam  was  also  tested  in  the  housing  on  the  duck  to  compare  how 
closely  it  resembles  a  clamped  boundary  condition.  The  FRF  and  coherence  for  the  in¬ 
flight  configuration  are  shown  below  in  Figure  32.  While  the  coherence  was  still 
excellent,  the  FRF  looks  significantly  different.  The  FRF  of  the  perfectly  clamped  beam 
looked  a  lot  cleaner  than  the  beam  in  the  duck  housing.  The  peaks  of  the  magnitude  plots 
have  been  reduced  significantly  signifying  that  the  housing  in  the  test  article  introduced 
damping  into  the  system.  The  beam  was  allowed  room  to  wiggle  and  energy  in  the 
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system  was  lost  in  this  manner.  The  looseness  of  the  clamp  also  reduced  the  frequency  of 
the  first  three  modes.  By  analyzing  the  eigenvectors  and  looking  for  the  first  three  beam 
bending  modes  the  first  two  frequencies  of  the  beam  were  found  and  shown  in  Table  4. 
The  third  mode  was  so  damped  on  the  FRF  that  the  actual  frequency  of  occurrence  could 
not  be  determined.  The  first  two  modes  dropped  by  23%  and  9%  respectively.  This  drop 
in  the  modes  signifies  that  the  duck  housing  was  more  compliant  than  an  ideal  clamp. 


Figure  32.  Comparison  of  the  measured  frequency  response  functions  of  the  beam  in  the 

test  article  and  clamped  configurations. 
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Table  4.  Comparison  of  measured  modes  for  the  leading  edge  beam  in  a  clamped 
configuration  and  in  the  test  article.  Frequencies  are  in  Hz. 


Beam  Configuration 

Mode  1 

Mode  2 

Mode  3 

Modulus  of 
Elasticity  (GPa) 

Clamped 

54.38 

336.9 

960 

- 

Abaqus  Model 

55.01 

345.0 

965.9 

88.1 

MATLAB  Clamped  Model  (untuned) 

55.07 

345.1 

966.3 

88.1 

MATLAB  Clamped  Model  (tuned) 

54.38 

340.8 

954.2 

85.9 

Experiment  in  Test  Article 

41.88 

306.3 

- 

- 

MATLAB  Loose  Clamp 

41.88 

279.6 

809.9 

85.9 

To  account  for  this  “wiggle  room”  in  the  physical  article,  another  optimization 
routine  was  developed  to  weaken  the  joint.  The  clamped  position  was  moved  toward  the 
duck  body  or  fuselage,  effectively  lengthening  the  free  end  of  the  beam  ever  so  slightly. 
An  element  was  created  equal  to  the  length  of  the  distance  the  pivot  point  was  moved. 
This  element  was  weakened  by  changing  the  length,  modulus  of  elasticity,  and  the  radius 
of  the  element.  This  element  provided  a  weak  link  that  attached  the  rest  of  the  bar  to  the 
clamp  in  order  to  mimic  the  actual  housing  the  beam  sat  in.  The  results  of  the 
optimization  are  shown  in  Table  5.  Since  the  original  length  of  the  beam  from  the 
clamped  end  was  135.8,  this  optimization  method  only  lengthened  the  beam  by  3.8% 
which  had  a  small  but  negligible  effect  on  the  overall  simulation.  The  advantage  of  this 
method  was  its  ease  of  implementation  inside  of  Abaqus. 
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Table  5.  Comparison  of  optimized  FE  model  versus  experimental  results  of  beam  in 

loose  housing. 


Mode  1 

Mode  2 

Mode  3 

In  Test  Article 

41.88 

306.3 

- 

MATLAB  Loose  Clamp 

41.88 

279.6 

809.9 

Optimization  Parameters 

Length 

(mm) 

Modulus 

(GPa) 

Radius 

(mm) 

Weakened  Beam  Element 

5.14 

27.1 

.453 

Normal  Beam  Element 

- 

85.9 

.487 

Modal  testing  was  also  perfonned  on  the  wing  to  assess  the  change  in  the  FRF  by 
adding  the  wing  to  the  beam.  Figure  33  shows  the  FRF  of  the  entire  wing — that  is  all  the 
FRFs  from  each  test  point  have  been  averaged  to  create  this  one  plot.  This  graph  gives 
the  overall  effect  of  the  membrane  on  the  beam.  As  expected  the  membrane  provided 
damping  to  the  beam  so  that  the  first  two  modes  occured  at  1 1.6  and  36.9  Hz.  The  beam 
in  a  clamped  configuration  had  the  first  mode  out  at  54  Hz  and  now  the  entire  wing  has 
two  modes  below  the  first  beam-bending  mode.  The  overall  FRF  presents  no  way  to 
distinguish  if  these  two  modes  are  membrane  or  beam  driven  modes.  As  such  the  FRFs 
of  six  distinct  points  on  the  wing  were  graphed  against  each  other  in  Figure  34.  These  six 
locations  were  the  beam  tip,  batten,  tape  in  the  middle  of  the  wing,  a  membrane  location 
in  the  middle  of  the  wing,  a  membrane  location  near  the  root  of  the  beam,  and  a 
membrane  location  near  the  aft  edge  of  the  wing.  The  exact  locations  of  these  points  are 
shown  pictorially  in  Figure  35. 
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Figure  33.  Overal  FRF  of  wing  in  clamped  configuration  and  in  the  test  article. 
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Figure  34.  FRF  and  coherence  for  wing  in  clamped  configuration  at  six  different  points. 
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Figure  35.  Locations  of  selected  FRF  plots  numbered  to  correspond  with  order  of 

appearance  in  the  legend. 


The  strongest  FRFs  at  the  first  mode  were  the  membrane,  edge  of  the  membrane,  and 
batten.  This  suggested  that  the  first  mode  of  the  wing  was  actually  driven  by  the 
membrane  and  not  the  beam.  Contrary,  at  the  second  mode,  the  beam  FRF  was  the 
highest  of  all  the  parts  so  this  mode  was  driven  by  the  beam  dynamics.  Without  the 
membrane,  the  clamped  beam  bending  mode  occurred  at  54.4  Hz  but  with  the  addition  of 
the  membrane  this  mode  has  been  reduced  by  32%  to  36.9  Hz.  Also,  while  the  coherence 
of  the  beam  was  generally  near  one  over  the  range  of  tested  frequencies,  the  coherence  of 
the  membrane  FRFs  dropped  below  acceptable  values  over  certain  frequency  ranges. 

The  coherence  was  acceptable  around  the  frequencies  of  the  first  two  modes  but  the 
coherence  of  the  membrane  dropped  below  0.8  for  several  frequency  bands  especially  in 
the  higher  frequencies.  Since  coherence  is  a  measure  of  the  linearity  between  the  input 
and  output,  either  the  membrane  did  not  behave  in  a  linear  manner  at  certain  frequencies 
or  the  data  recorded  at  some  points  were  poor.  It  was  probably  a  combination  of  both 
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explanations.  The  coherence  was  better  at  some  locations  versus  others  but  the  wing  was 
painted  to  give  a  good  return  to  the  laser  and  even  the  points  with  a  worse  coherence  still 
had  frequency  ranges  with  very  good  coherence.  Perhaps  the  membrane  did  not  behave 
as  linear  as  the  stiffer  materials  like  the  beam.  While  the  change  in  coherence  further 
shows  the  difficulties  in  modeling  membranes,  the  coherence  was  good  around  the 
frequencies  of  the  first  two  modes.  To  validate  the  effect  of  the  membrane  and  beam  on 
the  first  two  modes  the  mode  shapes  of  the  first  two  wing  modes  were  investigate  further. 

The  first  two  mode  shapes  of  the  wing  measured  in  the  lab  are  displayed  in  Figure  36. 
The  first  mode  at  1 1.6  Hz  was  characterized  by  high  membrane  deflections  especially  in 
the  region  highlighted  by  the  red  in  the  first  picture.  This  section  of  the  membrane  was 
farthest  from  the  clamped  boundary  condition  at  the  root  and  the  attachment  at  the 
leading  edge.  This  membrane  deflection  drove  the  first  mode  and  caused  the  beam  to 
deflect  but  the  beam  remained  mostly  rigid.  More  beam  deflection  and  bending  was  seen 
in  the  second  mode.  Here  the  maximum  velocity  occurred  out  at  the  tip  of  the  beam  and 
membrane.  The  section  of  the  membrane  previously  characterized  by  high  deflections  in 
the  first  mode  deflected  a  minute  amount  in  the  second  mode.  These  findings  confirmed 
that  the  membrane  drove  the  first  mode  and  the  beam  drove  the  second  mode.  The  beam 
shape  resembled  a  first  bending  mode.  Interestingly  enough,  the  duck  flaps  at  15  Hz  in 
air  which  is  closer  to  the  membrane  driven  mode  rather  than  the  first  beam  driven  mode. 
Previous  research  studies  have  shown  that  it  is  most  efficient  to  flap  near  the  resonance 
frequency  [7]  however  most  studies  have  been  performed  on  smaller  flexible  wings  or 
larger  more  rigid  wings  where  the  first  mode  is  a  beam  bending  mode.  Since  no  studies 
have  involved  the  use  of  highly  flexible  membrane,  either  the  flapping  frequency  of 
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highly  flexible  membrane  is  closer  to  the  first  bending  mode  or  this  design  was  simply 
not  optimized  to  flap  at  the  most  efficient  flapping  frequency.  The  latter  reason  was 
probably  true  because  this  COTS  SUAS  exhibits  poor  flying  qualities,  yet  there  could 
exist  a  correlation  between  this  membrane  mode  and  flapping  efficiency  but  without 
more  experimental  data  it  is  impossible  to  know  if  a  first  membrane  mode  is  more 
common  to  flexible  wings  or  simply  unique  to  this  test  article. 

Since  the  flapping  frequency  occurred  near  the  first  membrane  mode,  it  would  be 
expected  that  the  flapping  shape  would  look  similar  to  the  first  mode.  Air  and  vacuum 
testing  described  below  indeed  shows  that  the  wing  at  the  top  and  bottom  of  the  stroke 
closely  resembles  the  top  and  bottom  pictures  in  left  column  of  Figure  36.  However,  this 
laser  vibrometer  test  was  performed  with  the  wing  in  a  clamped  configuration.  Beam 
modal  testing  showed  there  was  a  difference  between  the  clamped  configuration  and  in 
the  test  article.  To  study  the  effect  of  the  boundary  condition  on  the  wing  FRF,  the  wing 
was  tested  attached  to  the  test  article — the  same  configuration  as  it  would  experience  in 
flight. 
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Figure  36.  Mode  shape  for  clamped  wing  at  1 1.6  Hz  (left)  and  36.9  Hz  (right) 
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The  same  six  points  graphed  in  Figure  34  above  are  compared  to  their  equivalent 
points  on  the  wing  attached  to  the  test  article  in  Figure  37.  The  phase  and  coherence  of 
these  points  can  be  found  in  the  Appendix.  As  with  the  beam,  the  FRFs  from  the  wing  on 
the  duck  were  less  clean  than  their  clamped  counterparts.  However,  looking  at  the  test 
point  on  the  beam,  the  max  peak  was  dampened  and  reduced  in  frequency  but  the  overall 
magnitude  in  the  system  was  relatively  similar.  The  first  mode  stayed  relatively 
unaffected  for  most  parts  of  the  wing.  Beyond  this  first  mode  the  FRF  was  altered  in 
different  ways  depending  on  the  part.  The  loose  clamp  introduced  damping  into  the 
system  but  since  it  left  the  first  mode  relatively  unchanged  the  mode  shape  was  expected 
to  be  similar  to  the  clamped  testing.  The  loose  nature  of  the  test  article  clamp  shows  that 
it  was  designed  to  operate  in  the  15  Hz  region.  Were  the  test  article  to  be  driven  at  a 
flapping  frequency  of  50  Hz  or  higher  the  wing  would  most  likely  interact  with  the 
housing  in  a  negative  way.  Also  the  housing  in  the  duck  has  seemed  to  significantly  alter 
the  FRFs  of  the  membrane  more  than  the  beam. 
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Figure  37.  Comparison  between  the  clamped  and  in  flight  configurations  for  six  different 

points  on  the  wing. 
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The  beam  FRF  changed  slightly  because,  although  the  boundary  condition  changed, 
the  beam  was  still  a  stiff  material.  On  the  other  hand,  the  other  boundary  condition  that 
changed  from  the  clamped  configuration  to  that  in  the  test  article  was  the  connector  piece 
on  the  back  of  the  membrane.  In  the  clamped  configuration,  the  connector  was  clamped 
between  the  same  two  plates  the  beam  was  clamped  between  so  that  there  was  no  slack  in 
the  membrane.  When  the  wing  was  attached  to  the  duck,  the  connector  was  angled 
slightly  below  the  beam  attachment  point.  This  created  a  slight  camber  in  the  wing 
membrane  as  seen  in  Figure  38.  This  camber  put  the  wing  in  a  different  state  of  tension 
compared  to  the  clamped  configuration  and  showed  that  the  membrane  was  depended  on 
the  initial  tension  applied  to  the  membrane.  This  initial  tension  made  the  membrane 
difficult  to  model  in  Abaqus. 


Figure  38.  Trace  of  the  camber  of  the  outside  of  the  wing. 

The  entire  duck  article  was  also  subjected  to  laser  vibrometer  testing.  This  testing 
gave  greater  insight  into  the  overall  modes  of  the  entire  system.  The  overall  test  article 
FRF  compared  to  just  the  wing  in  the  test  article  can  be  found  in  Appendix.  The  main 
difference  between  the  clamped  wing  and  the  entire  test  article  was  the  addition  of 
another  mode  at  15.6  Hz.  This  mode  was  cause  by  the  tail  and  characterized  by  a 
flapping  motion  seen  in  Figure  40.  The  first  mode  in  Figure  39  was  the  same  as  the 
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membrane  mode  described  above  and  likewise  the  mode  shape  displayed  in  Figure  41 
was  the  same  as  the  second  mode  in  the  clamped  wing  configuration.  A  summary  of  first 
three  wing  mode  frequencies  for  the  three  configurations  can  be  found  in  Table  6. 


Figure  39.  Mode  shape  at  10.63  Hz. 
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Figure  40.  Mode  shape  at  15.63  Hz. 
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Figure  41.  Mode  shape  at  33.28  Hz. 


Table  6.  Summary  of  laser  vibrometer  testing  for  the  wing. 


Frequencies  (Hz) 

Mode  1 

Mode  2 

Mode  3 

Wing  Clamped 

11.56 

36.88 

- 

Wing  In  Duck 

10.94 

31.25 

- 

Entire  Duck 

10.63 

15.63 

33.28 

4.5  Vacuum  Testing 

Model  tuning  gave  some  validation  that  the  FE  model  matched  the  test  article  but  this 
did  not  guarantee  that  the  FE  analysis  would  produce  the  same  dynamic  results  as  the 
wing  flapping  in  vacuum.  The  vacuum  test  setup  described  above  was  the  main  crux  and 
comparison  tool  used  to  validate  the  dynamic  simulations.  When  the  duck  was  hooked 
up  to  the  DC  power  supply  located  outside  the  vacuum  chamber  the  flapping  frequency 
of  the  duck  would  drop  from  15  to  12.5  Hz  even  though  the  manufacturer  specified  3.7 
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mV  were  supplied  to  the  duck.  This  small  change  in  frequency  did  not  matter  in  the  end 
because  ultimately  the  FE  model  could  be  modeled  to  reflect  any  excitation.  Vacuum 
testing  took  place  from  763  Torr  all  the  way  down  to  1.33  Torr.  The  low  pressure  of  the 
vacuum  chamber  provided  an  environment  that  proved  to  be  an  excellent  comparison 
against  the  FE  model  as  opposed  to  other  vacuum  tests  described  in  chapter  two  that 
either  utilized  a  vacuum  chamber  with  poor  view  of  the  test  article  or  did  not  have  access 
to  a  true  vacuum  chamber  and  had  to  replicate  the  effects  using  other  methods  such  as 
replacing  the  air  with  helium. 

The  flapping  frequency  at  a  range  of  pressures  was  evaluated  to  determine  the  effect 
of  the  vacuum  chamber  on  the  duck.  Figure  42  below  shows  that  as  pressure  decreases 
the  flapping  frequency  increases  in  an  exponential  manner.  An  exponential  curve, 
y  =  4. 5  9  A0'0053*  + 12.4  ,  fit  the  data  the  best  and  shows  that  at  the  lower  pressures  the 
flapping  frequency  was  more  sensitive  to  changes  in  pressure.  The  trend  shows  the 
importance  of  being  able  to  test  the  duck  in  a  proper  vacuum  chamber.  As  the  air  in  the 
chamber  was  removed  the  flapping  frequency  increased  for  a  constant  excitation  voltage 
because  less  damping  was  experienced  by  the  wing.  Since  the  wing  was  more  sensitive 
in  the  lower  pressure  region,  a  difference  of  30  Torr  corresponded  to  a  flapping  frequency 
reduction  of  approximately  1  Hz.  Therefore,  in  order  to  get  a  proper  comparison  between 
vacuum  data  and  a  FE  model,  a  relatively  small  change  in  air  pressure  could  cause 
significantly  different  behavior  in  the  wing  membrane. 
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Figure  42.  As  the  pressure  decreases  the  flapping  frequency  of  the  test  article  increases 

as  voltage  is  held  constant. 

One  area  of  interest  in  the  vacuum  test  was  the  qualitative  difference  between  the 
wing  behavior  in  air  versus  vacuum.  At  a  casual  glance  the  behavior  of  the  flapping 
motion  did  not  appear  significantly  different:  the  flapping  motion  stayed  constant  and  no 
irregularities  were  noticed.  Careful  inspection  of  the  high  speed  footage  revealed  the  true 
differences  between  was  responses  at  the  different  pressures.  Figure  43  and  Figure  44 
below  highlight  the  difference  between  the  flapping  strokes  at  atmospheric  pressure  (760 
Torr)  and  the  lowest  pressure  test  (1.33  Torr).  The  membrane  deflections  in  the  vacuum 
where  much  more  pronounced  than  in  air  as  was  expected.  Because  of  the  lack  of 
damping  in  the  absence  of  air  and  the  presence  of  larger  membrane  deflections,  the 
leading  edge  beam  deflected  farther  at  the  top  and  bottom  of  the  stroke  and  the  curvature 
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of  the  beam  was  more  pronounced.  This  extra  deflection  was  largely  due  to  the 
membrane.  In  air  the  leading  edge  of  the  wing  pulled  the  membrane  through  the  stroke 
and  at  the  end  of  each  stroke  the  momentum  of  the  membrane  caught  up  to  the  motion  of 
the  beam  so  that  the  beam  and  aft  part  of  the  membrane  arrived  at  the  bottom  of  the 
stroke  at  roughly  the  same  time.  In  vacuum,  however,  the  momentum  of  the  membrane 
was  not  slowed  by  aerodynamic  forces  and  as  a  result  when  the  beam  came  to  the  bottom 
of  the  stroke  the  membrane  still  had  extra  momentum  and  pulled  the  beam  even  farther  at 
the  top  and  bottom  of  the  stroke.  This  observed  membrane  behavior  was  a  key  feature  in 
the  qualitative  assessment  of  the  FE  model. 
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Figure  43.  High  speed  footage  captured  during  vacuum  testing.  The  two 
pictures  on  the  left  were  taken  at  atmospheric  conditions  (760  Torr).  The  two 
pictures  on  the  right  were  taken  at  1.33  Torr.  The  top  two  pictures  were  at  the 
top  of  the  stroke  while  the  bottom  two  pictures  were  the  bottom  of  the  stroke. 


87 


Figure  44.  Same  four  images  as  above  figure  overlaid  on  top  of  each  other.  The  beam 
and  membrane  deflections  of  the  near  vacuum  pressure  deflect  more  than  the  atmospheric 

tests. 


4.6  FE  Model 

Modeling  the  wing  in  Abaqus  proved  to  be  a  difficult  task.  While  FE  solvers  can 
excel  at  handling  stiff  lightly  damped  structures  the  flexible  nature  of  the  membrane 
brought  many  computational  hardships  with  it.  The  high  flexibility  of  the  membrane 
caused  the  FE  solver  to  predict  modes  approximately  every  5  Hz.  With  this  many  modes 
it  made  the  task  of  model  tuning  nearly  impossible  because  many  more  modes  were 
predicted  and  none  of  the  predicted  modes  matched  the  experimental  mode  shapes.  The 
response  of  the  wing  membrane  was  also  sensitive  to  the  input  and  material  properties. 
Because  of  all  these  reason,  validating  the  membrane  was  a  very  intricate  procedure  and 
tip  deflection  and  a  qualitative  assessment  of  the  deformations  were  used  as  primary 
methods  of  assessing  the  analysis. 

Many  models  were  run  in  Abaqus  to  find  the  parameters  that  produce  the  best 
dynamic  explicit  results.  For  modal  analysis  shell  elements  were  used  for  all  wing  parts 
(excluding  the  beam)  as  the  membrane  elements  proved  to  be  problematic  for  an 
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eigenvalue  analysis.  For  the  dynamic  simulations  several  iterations  were  run  with 
different  pennutations  of  shell  and  membrane  elements.  As  more  parts  were  converted  to 
shell  elements  the  computational  cost  rose.  Eventually,  the  final  configuration  was  found 
with  batten  and  connector  elements  as  shells  and  tape  and  membrane  as  membrane 
elements  as  detailed  in  Table  3.  This  configuration  was  the  happy  medium  between 
accuracy  and  computational  cost.  When  shell  elements  were  used  for  the  membrane 
portion  of  the  wing,  the  computational  cost  rose  significantly  and  the  membrane  behaved 
artificially  stiff  during  the  simulation.  All  simulations  were  run  using  Abaqus  6.1 1  on 
Linux  operating  system.  It  was  found  that  changing  the  damping  values  had  negligible 
effects  on  the  simulation  unless  they  increased  by  several  orders  of  magnitude. 

As  discussed  earlier,  the  wing  was  flapped  in  Abaqus  by  forcing  the  rotation  on  the 
beam.  This  forced  boundary  condition  gave  the  advantage  that  the  rotation  could  be 
explicitly  controlled.  The  downside  to  this  method  of  rotation  was  that  it  did  not 
accurately  capture  how  the  wing  was  actuated.  During  the  vacuum  testing  it  was  shown 
that  at  a  constant  voltage  the  flapping  frequency  increased  as  the  pressure  in  the  chamber 
decreased.  With  virtually  no  aerodynamic  forces  acting  on  the  wing  there  was  less 
resistance  and  the  flapping  frequency  increased.  By  forcing  the  rotation  at  the  root,  the 
FE  solver  would  not  be  able  to  repeat  this  behavior.  In  order  to  replicate  the  vacuum 
experiments  the  flapping  frequency  would  have  to  be  increased  to  17  Hz  from  the  15  Hz 
that  had  been  prescribed  in  all  previous  work  done  on  this  project.  Figure  45  shows  an 
early  comparison  of  the  full  wing  versus  the  Abaqus  explicit  bar  only  model  at  15  Hz 
flapping  frequency. 
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Figure  45.  Early  comparison  of  the  tip  deflection  versus  time  of  the  Abaqus  bar  only 
model  to  the  full  wing  model  flapping  at  15  Hz. 

One  of  the  most  apparent  differences  is  the  damping  of  the  transients.  The  bar  only 
model  exhibited  high  frequency  modes  on  top  of  the  15  Hz  flapping  frequency  as  well  as 
a  very  low  transient  increasing  in  amplitude  near  the  end  of  the  0.5  second  time  window. 
The  addition  of  the  wing  membrane  damped  out  these  high  frequency  transients.  A  low 
frequency  transient  was  still  observed  in  the  wing  model  but  the  amplitude  of  such  was 
negligible.  The  average  amplitudes  between  the  two  models  where  similar  but  it  was 
expected  that  the  wing  model  would  have  greater  deformations  due  to  the  addition  of  the 
membrane  alone.  In  the  vacuum  testing,  the  inertia  of  the  membrane  was  observed 
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creating  greater  deflections  during  stroke  reversal  but  this  was  not  seen  in  the  FE  element 
model.  To  further  compare  the  wing  model  to  the  vacuum  testing,  the  tip  deflection  was 
compared  to  the  vacuum  test  as  shown  in  Figure  46. 
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Figure  46.  Comparison  of  tip  deflection  of  the  FE  model  versus  vacuum  data  conformed 

to  15  Hz. 


It  is  important  to  note  that  the  time  vector  of  the  vacuum  data  has  been  conformed  to 


15  Hz  by  multiplying  it  by  ^  to  compare  it  against  the  FE  model.  The  actual  flapping 


frequency  of  the  vacuum  data  was  at  17  Hz.  For  now,  it  is  apparent  from  the  above 
figure  that  the  FE  solver  under  predicted  values  observed  in  the  vacuum  chamber.  The 
actual  tip  deflection  was  130  mm  but  the  FE  analysis  only  predicted  approximately  a  95 
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mm  deflection — a  difference  of  27%.  This  would  not  be  an  acceptable  difference  for 
modeling  purposes.  To  rectify  this  difference,  the  forced  flapping  rotation  in  the  FE 
model  was  reinvestigated.  Previously  models  of  the  wing  used  the  four  bar  linkage 
Simulink  simulation  which  outputted  a  stroke  amplitude  of  30°.  Physical  measurements 
of  the  test  article  showed  that  the  actual  stroke  amplitude  was  36.4°.  This  six  degree 
difference  at  the  root  could  account  for  a  large  difference  at  the  tip.  The  simulation  was 
rerun  with  this  new  stroke  amplitude  and  the  results  are  shown  below  in  Figure  47. 
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Figure  47.  Tip  deflection  comparison  with  FE  model  at  a  36.4  degree  stroke  amplitude 

with  vacuum  data  confonned  to  15  Hz. 


The  six  degree  increase  in  stroke  amplitude  has  recovered  approximately  12%  of  the 
amplitude  but  with  a  tip  deflection  amplitude  of  1 10  mm  the  FE  analysis  was  still  off  by 
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15%.  In  all  these  comparisons  the  time  vector  of  the  vacuum  data  has  been  conformed  to 
15  Hz,  so  the  forced  rotation  in  the  FE  model  was  increased  to  17  Hz  to  match  the  actual 
frequency  seen  during  the  vacuum  tests.  With  this  change  the  new  tip  deflection  is 
shown  below  in  Figure  48.  The  average  amplitude  increased  to  120  mm  corresponding  to 
a  8%  difference.  This  model  matched  the  vacuum  data  much  better  but  it  did  introduce 
more  transients  in  the  beginning.  The  tip  deflection  just  represents  one  physical  point  in 
the  entire  simulation  and  does  not  convey  how  the  membrane  acts. 

The  membrane  deflections  will  be  compared  qualitatively  using  the  high-speed 
images  from  vacuum  testing.  Figure  49  shows  a  comparison  of  the  down  and  upstroke  of 
between  the  FE  results  and  experimental  results.  The  most  obvious  difference  between 
the  two  on  the  down  stroke  is  the  membrane  deflection.  The  membrane  in  the 
experimental  testing  deflected  far  more  than  the  FE  model  at  both  the  top  and  the  bottom 
of  the  stroke.  Not  only  was  the  membrane  deflection  smaller  but  the  membrane 
dynamics  behaved  a  little  differently  than  the  experimental  results.  In  the  vacuum  testing 
results  the  membrane  can  be  seen  to  be  lagging  behind  the  leading  edge  beam  until  the 
wing  rotates  just  below  horizontal.  The  FE  program  predicted  the  membrane  leading  the 
beam  well  before  crossing  this  threshold  and  in  fact  the  membrane  led  more  than  the 
beam  during  the  down  stroke. 

The  upstroke  shows  similar  trends.  The  FE  membrane  deflections  were  less  than 
observed  except  during  the  end  of  the  stroke.  At  the  top  of  the  stroke  the  membrane 
deflections  looked  very  similar  and  the  point  in  the  stroke  where  the  membrane  crossed 
the  leading  edge  beam  was  similar  to  the  observed  results.  The  upstroke  on  the  test 
article  happened  faster  than  that  of  the  FE  predictions.  This  was  evident  by  the  uneven 
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spacing  of  the  wing  positions  of  the  vacuum  data.  In  fact,  the  FE  element  stroke  looked 
symmetrical — the  down  stroke  and  upstroke  were  almost  identical.  This  was  not  the  case 
with  the  vacuum  results.  The  membrane  deflections  were  different  in  the  upstroke  and 
down  stroke  and  the  moment  the  membrane  crossed  the  beam  occurred  at  different 
phases  during  the  up  and  down  stroke.  The  difference  in  the  flapping  motions  is  apparent 
in  Figure  50. 
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Figure  48.  Tip  deflection  comparison  with  FE  model  run  at  17  Hz. 
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Figure  49.  Comparison  of  FE  simulation  versus  vacuum  testing.  Top  row  shows  the 
down  stroke,  bottom  row  shows  the  upstroke. 


Figure  49  shows  the  tip  deflection  from  eight  flap  cycles  plotted  on  top  of  one 
another.  The  FE  data  was  taken  from  0.5-1  seconds  after  the  run  began  to  mitigate  the 
effect  of  transients.  The  vacuum  data  represents  the  average  steady  state  values  measured 
in  the  lab.  The  blue  line  represents  the  shape  of  the  input  signal  going  into  the  FE 
simulation.  The  amplitude  and  bias  of  the  input  signal  have  been  altered  to  match  the 
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scale  of  the  graph  so  the  y-axis  values  corresponding  to  the  blue  line  are  not  meaningful, 
but  what  is  meaningful  is  the  shape  of  the  input.  The  input  signal  was  not  a  pure 
sinusoid — close  but  not  a  sinusoid. 

The  FE  data  was  relatively  repeatable  and  would  converge  toward  a  steady  state 
solution  if  run  for  a  long  enough  simulation  time.  The  apparent  difference  in  the  same  17 
Hz  stroke  comes  from  the  upstroke.  The  Abaqus  results  overall  look  symmetric  and 
sinusoidal.  The  input  signal  is  shown  for  a  frame  of  reference  and  not  to  suggest  that  the 
tip  deflection  should  follow  the  input  signal  because  it  should  not.  The  tip  of  a  pure  rigid 
body  should  display  the  same  signal  shape  as  the  input  but  because  of  the  deflections  of 
the  beam  the  actual  dip  path  will  deviate.  The  important  takeaway  is  that  the  input  to  the 
FE  element  program  was  not  purely  sinusoidal  but  the  output  at  the  tip  was  more 
sinusoidal.  This  differs  greatly  from  the  tip  deflection  of  the  vacuum  data.  The  upstroke 
occurs  faster  than  the  down  stroke  causing  a  skewed  sinusoid.  This  feature  was  not 
observed  in  the  FE  simulation  and  was  probably  mostly  likely  due  to  an 
oversimplification  in  the  input  mechanisms  into  the  FE  simulation. 
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Figure  50.  Tip  deflection  from  eight  flap  cycles.  The  FE  data  comes  from  after  0.5 
seconds  into  the  run.  The  vacuum  data  is  an  average  of  the  steady  state  values. 


The  FE  model  with  the  weakened  joint  described  earlier  was  run  in  Abaqus.  While 
the  method  worked  for  tuning  the  model,  it  failed  to  produce  a  workable  model  for  the 
dynamic  simulation.  The  weakened  joint  was  simply  too  weak  to  withstand  the  forces 
imparted  by  the  rotation  of  the  wing.  As  a  result,  the  joint  defonned  excessively 
producing  unusable  data.  Tuning  the  model  with  the  clamped  data  proved  to  be  the  more 
accurate  and  reliable  method.  Had  the  weakened  joint  worked  properly,  it  should  have 
captured  a  small  portion  of  the  behavior  of  the  physical  model. 
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V. 


Conclusions  and  Future  Work 


5.1  Overview 

Modeling  the  flapping  wing  motion  remains  an  area  of  ongoing  study  because  of  the 
importance  to  aircraft  designers.  With  an  understanding  of  the  inertial  and  aerodynamic 
forces  aircraft  designers  could  optimize  SUAS  designs  to  meet  the  rising  demand  for 
smaller  operational  vehicles.  Most  research  has  been  aimed  at  figuring  out  the  low 
Reynolds  number  and  unsteady  aerodynamic  effects  that  characterized  flapping  wing 
flight.  Less  research  has  focused  on  validating  structural  models  of  flapping  wings.  As 
wings  become  more  flexible,  the  structural  response  becomes  more  nonlinear  and  more 
effort  must  be  focused  on  structural  validation. 

A  FE  model  of  a  flapping  wing  from  a  COTS  SUAS  was  made  and  compared  against 
the  actual  test  article.  Model  tuning  brought  the  first  beam  bending  of  the  leading  edge 
beam  of  the  wing  to  within  1%  of  the  test  article.  Wing  modal  tuning  was  not 
accomplished  due  to  the  myriad  of  modes  predicted  in  the  FE  element  solver.  The  FE 
solver  had  difficulties  with  the  large  flexible  membrane  area  and  none  of  the  predicted 
modes  matched  the  measured.  Wing  modal  testing  showed  that  the  addition  of  the 
membrane  to  the  beam  decreased  frequency  the  first  beam  bending  mode  by  32%.  The 
first  wing  mode  was  found  to  be  driven  by  the  membrane  and  not  the  beam.  This  mode 
was  closest  to  the  flapping  frequency  and  the  mode  shape  resembled  the  flapping  motion 
of  the  test  article.  The  tip  deflection  amplitude  was  predicted  to  within  8%  of  vacuum 
data.  Overall  predicted  membrane  deflections  were  lower  than  those  observed  in  the 
vacuum  chamber. 
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5.2  Significance  of  Research 

This  research  ventured  into  an  area  of  flapping  wing  studies  that  has  been  relatively 
untouched.  No  significant  research  has  delved  into  the  structural  aspect  of  a  billowing 
flapping  wing.  With  no  prior  examples  to  base  the  work  upon,  this  research  encountered 
many  unforeseen  complexities  involved  with  highly  flexible  structures.  The  complex 
nature  of  the  membrane  dynamics  required  precise  modeling  of  the  problem.  Everything 
from  the  dimensions,  material  properties,  and  input  signals  were  had  significant  influence 
on  the  results. 

While  the  modulus  of  elasticity  and  density  were  both  important  parameters 
controlling  how  the  membrane  deflected,  the  membrane  density  had  the  most  significant 
impact  on  the  membrane  kinematics.  A  small  change  to  the  membrane  density  changed 
the  transient  response  and  tip  deflection.  Increasing  the  density  increased  the  mass  of  the 
membrane  which  not  only  decreased  the  modes  at  which  the  membrane  was  excited  at 
but  the  added  mass  made  the  membrane  more  dominate  in  the  beam  membrane 
relationship  in  regards  to  flapping.  When  the  membrane  had  more  mass  it  could  cause 
greater  tip  deflections  because  it  had  more  influence  over  the  lighter  beam.  However,  the 
mass  and  stiffness  had  a  coupling  effect.  If  the  membrane  excitation  was  not  in  sync  with 
the  beam  then  a  heavier  membrane  could  cause  less  tip  deflection  because  the  membrane 
motion  would  oppose  the  beam  motion.  In  this  manner,  the  density  and  modulus  of 
elasticity  were  two  important  sensitive  parameters  toward  exacting  the  FE  response. 

Through  modal  testing,  the  effect  of  the  addition  of  a  membrane  to  a  leading  edge 
spar  was  seen.  Not  only  did  the  membrane  decrease  the  beam  natural  frequency  by  32% 
but  the  first  membrane  mode  was  closest  to  the  actual  flapping  frequency.  It  was  shown 
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that  modeling  a  loose  joint  by  weakening  the  material  properties  worked  for  modal 
analysis  but  proved  useless  for  dynamic  simulation. 

5.3  Future  Work 

While  the  general  flapping  nature  predicted  by  the  FE  program  looked  to  be 
fundamentally  the  same  motion  as  produced  by  the  test  article,  the  differences  in  flapping 
kinematics  could  be  improved  upon.  Reasons  for  different  membrane  kinematics  could 
have  resulted  from  the  approximations  made  in  the  FE  solver.  Membrane  elements 
where  used  for  the  Mylar  membrane  portion  instead  of  plates.  It  was  found  that  shell 
elements  worked  well  for  the  stiffer  parts  of  the  wing  but  if  used  for  the  highly  flexible 
membrane  then  the  shell  elements  made  it  unnaturally  stiff.  The  density  of  the  membrane 
could  be  increased  or  the  modulus  of  elasticity  decreased  in  an  attempt  to  get  more 
deflection  out  of  the  membrane  but  these  values  are  difficult  to  tune  because  the  dynamic 
simulations  can  be  computationally  intensive. 

Another  culprit  could  potentially  be  the  forced  rotation  at  the  root.  It  was  shown 
above  how  just  a  small  change  in  the  stroke  amplitude  or  the  flapping  frequency 
increased  the  tip  deflection.  The  previous  four  bar  Simulink  simulation  was  used  as  an 
input  for  every  test  case.  While  the  output  of  the  simulation  was  not  a  pure  sinusoidal 
signal  the  FE  flapping  motions  look  more  sinusoidal  than  the  test  article  data.  The  four 
bar  linkage  could  have  overly  simplified  the  input.  A  small  change  in  the  input  had 
significant  changes  in  the  tip  deflections  as  explained  in  Chapter  4.  The  same  was 
probably  true  of  the  membrane  kinematics.  Figuring  out  the  actual  input  to  the  wing 
warrants  further  investigation  by  taking  measured  data  rather  than  relying  on  simulations. 
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Another  possible  cause  of  the  difference  in  wing  kinematics  and  tip  deflection  could 
have  been  the  location  of  the  input.  For  all  the  models,  the  input  location  was  assumed  to 
be  where  the  membrane  met  the  leading  edge  beam  (the  same  location  as  previous 
studies).  A  closer  look  at  the  wing  in  the  housing  revealed  a  slight  gap  (less  than  2  mm) 
between  the  wing  attachment  point  and  the  pivot  point.  While  not  a  significant  distance, 
this  gap  would  cause  a  longer  rotation  arm  and  could  account  for  some  difference  in  tip 
deflection  as  well  as  the  wing  kinematics. 

The  low  stiffness  of  the  membrane  made  the  model  difficult  to  tune.  Focusing  solely 
on  modal  tuning  of  a  flexible  membrane  model  would  greatly  benefit  the  flapping  wing 
field.  Modal  testing  allowed  the  tuning  of  the  leading  edge  beam  but  the  membrane  was 
not  able  to  be  tuned.  The  FE  solver  predicted  too  many  mode  shapes — none  of  which 
matched  the  experimental  results.  Finding  a  way  to  tune  the  membrane  would  give  more 
confidence  in  the  dynamic  simulation.  The  tension  on  the  membrane  plays  a  large  part  in 
this  model  tuning.  To  what  degree  is  unknown  but  the  lack  of  tension  on  the  membrane 
in  the  FE  model  could  have  made  a  difference  in  the  modal  analysis. 

Initial  intentions  were  to  supplement  the  FE  model  with  aerodynamic  forces  in  an 
attempt  to  realize  an  aeroelastic  model  that  did  not  rely  on  CFD.  Abaqus  allows  for  a 
subroutine  to  define  a  force  field  that  changes  with  each  subroutine.  A  combination  of 
BET,  empirical  methods,  and  unsteady  models  could  be  implemented  into  the  user 
subroutine  to  add  aerodynamic  loading  to  the  structure.  While  an  initial  shell  of  the  code 
needed  to  implement  the  above  method  was  developed,  time  constraints  hindered  full 
development.  Future  efforts  could  easily  leverage  the  existing  framework  to  add  a 
complete  aerodynamic  model  to  the  system. 
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While  mimicking  flapping  wing  motions  in  a  FE  solver  is  possible,  it  is  not  a  trivial 
task.  The  wing  membrane  was  sensitive  to  small  changes  in  density,  input  conditions, 
and  modulus  of  elasticity.  Although  running  a  dynamic  simulation  was  not  as 
computationally  intensive  as  coupling  a  CFD  and  FE  solver  together,  the  task  took  long 
enough  to  make  running  different  permutations  difficult.  It  is  recommended  that  any 
future  work  on  this  project  leverage  the  use  of  the  network  of  clustered  processors.  This 
would  reduce  the  computation  time  and  greatly  speed  up  the  process. 

5.4  Summary 

The  foundation  for  modeling  a  highly  flexible  flapping  wing  was  established.  The 
flapping  motion  and  dynamic  simulation  can  be  complex  so  precise  modeling  of  the 
flapping  mechanism  and  wing  geometry  is  extremely  important.  Everything  from 
geometric  measurements,  material  properties,  boundary  conditions,  and  activation  plays 
an  essential  part  toward  predicting  the  exact  behavior.  The  flexible  membrane  proved  to 
be  the  most  difficult  part  to  model  in  the  wing.  While  incorporating  the  more  subtle 
aspects  of  the  problem  into  the  model  would  help  to  better  results,  membrane  modal 
analysis  and  dynamic  response  will  still  prove  difficult.  Future  efforts  should  focus  more 
on  correctly  modeling  the  membrane  portion.  More  advance  analysis  techniques  within 
Abaqus  should  be  investigate  because  the  standard  membrane  and  shell  elements  have 
difficulties  predicting  the  nonlinearity  of  the  wing  membrane.  Future  efforts  should 
focus  on  improving  membrane  modal  testing,  updating  the  established  model  to 
incorporate  the  more  subtle  aspects  of  the  problem,  and  add  aerodynamic  forces  into  the 
model  through  a  user  subroutine. 
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Figure  A-  1.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 

and  in  the  test  article  at  the  tip  of  the  beam. 
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Figure  A-  2.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 

and  in  the  test  article  at  the  batten. 
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Figure  A-  3.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 
and  in  the  test  article  at  the  middle  of  the  tape  area. 


Membran  middle 


-50 

-100 

200 

0 

-200 


50 


100 


150 


200 


250 


'll 


Is 


100 


150 


200 


250 


Figure  A-  4.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 
and  in  the  test  article  at  the  middle  of  the  membrane. 
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Figure  A-  5.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 
and  in  the  test  article  at  the  membrane  near  the  base  of  the  wing. 
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Figure  A-  6.  Frequency  response  function  comparison  of  wing  in  clamped  configuration 
and  in  the  test  article  at  the  edge  of  the  membrane. 
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Figure  A-  7.  Wing  averaged  frequency  response  function  comparison  between  clamped 

configuration  and  whole  test  article. 
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